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)
src_fd = in_defn.GetFieldDefn(fld_index)
fd = ogr.FieldDefn(src_fd.GetName(), src_fd.GetType())
fd.SetWidth(src_fd.GetWidth())
fd.SetPrecision(src_fd.GetPrecision())
shp_layer.CreateField(fd)
#############################################################################
# Apply spatial query.
in_layer.SetSpatialFilterRect(s_minx, s_miny, s_maxx, s_maxy)
#############################################################################
# Setup rect geometry to try intersect with.
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint_2D(s_minx, s_miny)
ring.AddPoint_2D(s_maxx, s_miny)
ring.AddPoint_2D(s_maxx, s_maxy)
ring.AddPoint_2D(s_minx, s_maxy)
ring.AddPoint_2D(s_minx, s_miny)
filt_geom = ogr.Geometry(ogr.wkbPolygon)
filt_geom.AddGeometry(ring)
#############################################################################
# Process all features in input layer.
in_feat = in_layer.GetNextFeature()
while in_feat is not None:
geom = in_feat.GetGeometryRef()
if not is_paleo:
feature.SetField(i, str(timestamp))
polygon = None
# Save feature
layer.CreateFeature(feature)
# Cleanup
feature = None
else:
for level in contour_levels:
contour_var = np.array(np.squeeze(ppt.permute(nc.variables[varname], var_order)), order="C")
contour_points = get_contours(contour_var, x, y, nc_projection, level)
# For each contour
polygon = ogr.Geometry(ogr.wkbPolygon)
for k in range(0, len(contour_points)):
geoLocations = contour_points[k]
ring = ogr.Geometry(ogr.wkbLinearRing)
# For each point,
for pointIndex, geoLocation in enumerate(geoLocations):
ring.AddPoint(geoLocation[0], geoLocation[1])
ring.CloseRings()
polygon.AddGeometry(ring)
# Create feature
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(polygon)
feature.SetFID(k)
i = feature.GetFieldIndex("level")
feature.SetField(i, level)
polygon = None
# Save feature
layer.CreateFeature(feature)
# Cleanup
def SetZ (poGeom, dfZ ):
if poGeom is None:
return
eGType = wkbFlatten(poGeom.GetGeometryType())
if eGType == ogr.wkbPoint:
poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ)
elif eGType == ogr.wkbLineString or \
eGType == ogr.wkbLinearRing:
for i in range(poGeom.GetPointCount()):
poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ)
elif eGType == ogr.wkbPolygon or \
eGType == ogr.wkbMultiPoint or \
eGType == ogr.wkbMultiLineString or \
eGType == ogr.wkbMultiPolygon or \
eGType == ogr.wkbGeometryCollection:
for i in range(poGeom.GetGeometryCount()):
SetZ(poGeom.GetGeometryRef(i), dfZ)
geoTransform = raster.GetGeoTransform()
sourceSRS = osr.SpatialReference()
sourceSRS.ImportFromWkt( raster.GetProjectionRef() )
targetSRS = osr.SpatialReference()
targetSRS.ImportFromEPSG(4326)
coordTrans = osr.CoordinateTransformation(sourceSRS,targetSRS)
features = []
# simple sliding detection window
y0 = 0
while h-y0 >= winY:
x0 = 0
while w-x0 >= winX:
# Create geometry
ring = ogr.Geometry(ogr.wkbLinearRing)
xc,yc = tfRasterToProj(x0,y0,geoTransform); ring.AddPoint( xc,yc )
xc,yc = tfRasterToProj(x0+winX,y0,geoTransform); ring.AddPoint( xc,yc )
xc,yc = tfRasterToProj(x0+winX,y0+winY,geoTransform); ring.AddPoint( xc,yc )
xc,yc = tfRasterToProj(x0,y0+winY,geoTransform); ring.AddPoint( xc,yc )
xc,yc = tfRasterToProj(x0,y0,geoTransform); ring.AddPoint( xc,yc )
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# Transform to target SRS
poly.Transform(coordTrans)
# Read data
data = raster.ReadAsArray(x0, y0, winX, winY).astype(np.float)
# Classify data. Now this depends on if there is one or many feature vectors being computed
# handle those cases accordingly, maybe a majority decision, maybe count labels, etc
for featureVector in computeFeatureVectors(data):
labels = classifier.predict(featureVector)
layer.CreateField(ogr.FieldDefn('ROW', ogr.OFTInteger))
layer.CreateField(ogr.FieldDefn('COL', ogr.OFTInteger))
# Strata field
layer.CreateField(ogr.FieldDefn('STRATUM', ogr.OFTInteger))
# Loop through samples adding to layer
for i, (stratum, col, row) in enumerate(zip(strata, cols, rows)):
# Feature
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField('ID', int(i))
feature.SetField('ROW', int(row))
feature.SetField('COL', int(col))
feature.SetField('STRATUM', int(stratum))
# Geometry
ring = ogr.Geometry(type=ogr.wkbLinearRing)
for corner in corners:
ring.AddPoint(
gt[0] + (col + corner[0]) * gt[1] + (row + corner[1]) * gt[2],
gt[3] + (col + corner[1]) * gt[4] + (row + corner[1]) * gt[5])
square = ogr.Geometry(type=ogr.wkbPolygon)
square.AddGeometry(ring)
feature.SetGeometry(square)
layer.CreateFeature(feature)
feature.Destroy()
sample_ds = None
elif EQUAL(args[iArg], "-clipsrcsql") and iArg < nArgc - 1:
pszClipSrcSQL = args[iArg + 1]
iArg = iArg + 1
elif EQUAL(args[iArg], "-clipsrclayer") and iArg < nArgc - 1:
pszClipSrcLayer = args[iArg + 1]
iArg = iArg + 1
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(
def discover_geom_name(ogr_type):
"""
:param ogr_type: ogr GetGeomType()
:return: string geometry type name
"""
return {ogr.wkbUnknown : "UNKNOWN",
ogr.wkbPoint : "POINT",
ogr.wkbLineString : "LINESTRING",
ogr.wkbPolygon : "POLYGON",
ogr.wkbMultiPoint : "MULTIPOINT",
ogr.wkbMultiLineString : "MULTILINESTRING",
ogr.wkbMultiPolygon : "MULTIPOLYGON",
ogr.wkbGeometryCollection : "GEOMETRYCOLLECTION",
ogr.wkbNone : "NONE",
ogr.wkbLinearRing : "LINEARRING"}.get(ogr_type)
a += 1
a = 0
b += 1
return verts
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")
elif EQUAL(args[iArg], "-s_srs") and iArg < nArgc - 1:
iArg = iArg + 1
pszSourceSRSDef = args[iArg]
elif EQUAL(args[iArg], "-a_srs") and iArg < nArgc - 1:
iArg = iArg + 1
pszOutputSRSDef = args[iArg]
elif EQUAL(args[iArg], "-t_srs") and iArg < nArgc - 1:
iArg = iArg + 1
pszOutputSRSDef = args[iArg]
bTransform = True
elif EQUAL(args[iArg], "-spat") and iArg + 4 < nArgc:
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]))
poSpatialFilter = ogr.Geometry(ogr.wkbPolygon)
poSpatialFilter.AddGeometry(oRing)
iArg = iArg + 4
elif EQUAL(args[iArg], "-where") and iArg < nArgc - 1:
iArg = iArg + 1
pszWHERE = args[+ +iArg]
elif EQUAL(args[iArg], "-select") and iArg < nArgc - 1: