Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def sasplanet_hlg2ogr(fname):
with open(fname) as f:
lines = f.readlines(4096)
if not lines[0].startswith('[HIGHLIGHTING]'):
return None
coords = [[], []]
for l in lines[2:]:
val = float(l.split('=')[1].replace(',','.'))
coords[1 if 'Lat' in l else 0].append(val)
points = zip(*coords)
ld('sasplanet_hlg2ogr', 'points', points)
ring = ogr.Geometry(ogr.wkbLinearRing)
for p in points:
ring.AddPoint(*p)
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring)
ds = ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' )
assert ds is not None, 'Unable to create datasource'
#~ src_srs = txt2srs('EPSG:4326')
src_srs = txt2srs('+proj=latlong +a=6378137 +b=6378137 +datum=WGS84 +nadgrids=@null +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs')
layer = ds.CreateLayer('sasplanet_hlg', srs=src_srs)
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetGeometry(polygon)
layer.CreateFeature(feature)
del feature
del polygon
print("Shifting vector -> {}".format(os.path.basename(outputVectorFile)))
outVector = outDriver.CreateDataSource(outputVectorFile)
outSrs = osr.SpatialReference(outProjection)
# create layer
outLayer = outVector.CreateLayer(os.path.basename(outputLayerName),
srs=outSrs, geom_type=ogr.wkbPolygon)
outFeatureDef = outLayer.GetLayerDefn()
# create rings from input rings by shifting points
for feature in inputFeatures:
# create the poly
outPoly = ogr.Geometry(ogr.wkbPolygon)
poly = feature.GetGeometryRef()
for ring_idx in range(poly.GetGeometryCount()):
ring = poly.GetGeometryRef(ring_idx)
# create the ring
outRing = ogr.Geometry(ogr.wkbLinearRing)
for i in range(0, ring.GetPointCount()):
pt = ring.GetPoint(i)
outRing.AddPoint(pt[0] + offsetGeo[0], pt[1] + offsetGeo[1])
outPoly.AddGeometry(outRing)
# create feature
outFeature = ogr.Feature(outFeatureDef)
outFeature.SetGeometry(outPoly)
outLayer.CreateFeature(outFeature)
def world_to_pixel_poly(rpc_dict, geometry):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
pixelRing = ogr.Geometry(ogr.wkbLinearRing)
geoRing = geometry.GetGeometryRef(0)
numPoints = geoRing.GetPointCount()
rpc = rpc_from_gdal_dict(rpc_dict)
for p in range(numPoints):
point = numpy.array(map(float, geoRing.GetPoint(p)))
pixel, line = rpc.project(point)
pixelRing.AddPoint(pixel, line)
pixelPoly = ogr.Geometry(ogr.wkbPolygon)
pixelPoly.AddGeometry(pixelRing)
return pixelPoly
if os.path.exists(path):
os.remove(path)
outDriver = ogr.GetDriverByName('ESRI Shapefile')
outDataSource = outDriver.CreateDataSource(path)
outLayer = outDataSource.CreateLayer(path, geom_type=ogr.wkbPolygon)
featureDefn = outLayer.GetLayerDefn()
self.vert_list = vert_list(xs,ys)
for vert_set in self.vert_list:
ring = ogr.Geometry(ogr.wkbLinearRing)
first_point = vert_set[0]
for point in vert_set:
ring.AddPoint(point[0], point[1])
ring.AddPoint(first_point[0], first_point[1])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(poly)
outLayer.CreateFeature(outFeature)
outFeature.Destroy()
if verts ==True:
import gisexport
gisexport.export_coords_to_shp(self.vert_list, "gridverts.shp")
angles[i-1] = segments_angle(points[i-2], points[i-1], points[i])
if i == p_len-2:
angles[0] = segments_angle(points[p_len-2], points[0], points[1])
else:
angles[i+2] = segments_angle(points[i+1], points[i+2], points[i+3])
p_len = p_len + 1
else:
angles[i] = 2*MIN_ANGLE #angle was sharp but segments were too long, skip it in next iteration
#if segments are smaller then threshold
# if smallest angle is on starting point
i = np.argmin(angles) #index of sharpest angle
min_a = np.amin(angles) #sharpest angle
#while angle is smaller then threshold
if geom_type == 0:
ring = ogr.Geometry(ogr.wkbLinearRing)
else:
ring = ogr.Geometry(ogr.wkbLineString)
for i in range (0,p_len):
ring.AddPoint(points[i,0], points[i,1])
if points[0,0]!=points[p_len-1,0] or points[0,1]!=points[p_len-1,1]:
ring.AddPoint(points[0,0], points[0,1])
return 1,ring
#output: gen_ring=1 indicates that some geometry is preserved after generalisation
#output: out_geom will receive all generalised rings in it
#input: True means that rings are allways closed
#input: first 0 means that ogrLinearRing is geometry type for creation
#input: out_geom is reference to geometry in which new generalised geometry of right type will be stored
#input: alg_type - generalisation (simplification+smoothing), simplification or smoothing
#input: second 0 means that this is multigeometry (polygon with multiple rings)
gen_ring,out_geom = Decide(ring,True,0,out_geom,alg_type,0) #perform generalisation
#there were some geometry preserved, store this in gen so output feature will be created
#some parts could be deleted due to too small area for given scale
if gen_ring == 1:
gen = 1
elif geom.GetGeometryName() == 'MULTILINESTRING':
out_geom = ogr.Geometry(ogr.wkbMultiLineString) #create output geometry of given type
for i in range(0, geom.GetGeometryCount()): #iterate over lines
line = geom.GetGeometryRef(i)
#check if it closed polyline, if so, generalize it as ring, which means
#that neither vertex is considered as fixed
#if line is open, starting and ending points are preserved
ps = line.GetPoint(0)
pe = line.GetPoint(line.GetPointCount()-1)
closed = False
if (ps[0] == pe[0]) and (ps[1] == pe[1]):
closed = True
#output: gen_line=1 indicates that some geometry is preserved after generalisation
#output: out_geom will receive all generalised lines in it
#input: closed indicates whether line is closed or not
#input: 1 means that ogrLineString is geometry type for creation
#input: out_geom is reference to geometry in which new generalised geometry of right type will be stored
def to_ogr_multi_type(self, geom, ogr_type):
"""
Convert ogr geometry to multi-type when supplied the ogr.type.
:param geom: The ogr geometry
:type geom: OGRGeometry
:param ogr_type: The ogr geometry type
:type ogr_type: String
:return: WkB geometry and the output layer geometry type.
:rtype: Tuple
"""
multi_geom = ogr.Geometry(ogr_type)
multi_geom.AddGeometry(geom)
geom_wkb = multi_geom.ExportToWkt()
geom_type = multi_geom.GetGeometryName()
return geom_wkb, geom_type
elif EQUAL(args[iArg], "-clipsrcwhere") and iArg < nArgc - 1:
pszClipSrcWhere = args[iArg + 1]
iArg = iArg + 1
elif EQUAL(args[iArg], "-clipdst") and iArg < nArgc - 1:
if IsNumber(args[iArg + 1]) and iArg < nArgc - 4:
oRing = ogr.Geometry(ogr.wkbLinearRing)
oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 2]))
oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 4]))
oRing.AddPoint_2D(float(args[iArg + 3]), float(args[iArg + 4]))
oRing.AddPoint_2D(float(args[iArg + 3]), float(args[iArg + 2]))
oRing.AddPoint_2D(float(args[iArg + 1]), float(args[iArg + 2]))
poClipDst = ogr.Geometry(ogr.wkbPolygon)
poClipDst.AddGeometry(oRing)
iArg = iArg + 4
elif (len(args[iArg + 1]) >= 7 and EQUAL(args[iArg + 1][0:7], "POLYGON") ) or \
(len(args[iArg + 1]) >= 12 and EQUAL(args[iArg + 1][0:12], "MULTIPOLYGON")):
poClipDst = ogr.CreateGeometryFromWkt(args[iArg + 1])
if poClipDst is None:
print(
"FAILURE: Invalid geometry. Must be a valid POLYGON or MULTIPOLYGON WKT\n")
return Usage()
iArg = iArg + 1
elif EQUAL(args[iArg + 1], "spat_extent"):
iArg = iArg + 1
def make_ogr_point(x,y):
return ogr.Geometry(wkt="POINT(%f %f)"%(x,y))
for i in range(1, p_len):
p = g.GetPoint(i)
#remove duplicate neighbour points
if points[j,0] != p[0] or points[j,1] != p[1]:
j = j+1
points[j,0] = p[0]
points[j,1] = p[1]
p_len = j+1 #final number of points in point list
#malformed closed geometry
if p_len < 4:
return 0,g
#triangle, nothing to generalize, return cleaned
if p_len == 4:
if geom_type == 0:
ring = ogr.Geometry(ogr.wkbLinearRing)
else:
ring = ogr.Geometry(ogr.wkbLineString)
for i in range (0,p_len):
ring.AddPoint(points[i,0], points[i,1])
if points[0,0]!=points[p_len-1,0] or points[0,1]!=points[p_len-1,1]:
ring.AddPoint(points[0,0], points[0,1])
return -1,ring
#calculate squared segments lengths, algorithm always change line where shorthest segment is found
segments = np.zeros(p_len-1)
for i in range(0, p_len-1):
segments[i] = squared_length(points[i],points[i+1])
i = np.argmin(segments) #index of shorthest segment
min_s = np.amin(segments) #squared length of shorthest segment
while min_s < SQR_TOL_LENGTH and p_len > 4: