def convertNodesCoordinates(graph, sourceEPSG, targetEPSG):
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(sourceEPSG)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(targetEPSG)
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
for node in graph.nodes():
x, y = graph.node[node]['longitude'], graph.node[node]['latitude']
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x,y)
point.Transform(coordTrans)
graph.node[node]['longitude'], graph.node[node]['latitude'] = point.GetX(), point.GetY()
return graph
python类Geometry()的实例源码
def convertCoordinates(coordinate, sourceEPSG, targetEPSG):
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(sourceEPSG)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(targetEPSG)
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
point = ogr.Geometry(ogr.wkbPoint)
lon, lat = coordinate
point.AddPoint(lon,lat)
point.Transform(coordTrans)
return point.GetX(), point.GetY()
def batchConvertCoordinates(sourceCoordinates, sourceEPSG, targetEPSG):
targetCoordinates = []
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(sourceEPSG)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(targetEPSG)
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
for lon, lat in sourceCoordinates:
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(lon,lat)
point.Transform(coordTrans)
targetCoordinates.append((point.GetX(), point.GetY()))
return targetCoordinates
def get_convex_hull(shape_name, destination_directory):
'''
This method will read all objects from a shape file and create a new one
with the convex hull of all the geometry points of the first.
'''
driver = ogr.GetDriverByName(str('ESRI Shapefile'))
shape = driver.Open(shape_name, 0)
layer = shape.GetLayer()
layer_name = layer.GetName()
spatial_reference = layer.GetSpatialRef()
prefix = get_basename(shape_name)
output_name = create_filename(destination_directory, '%s-hull.shp' % prefix)
geometries = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in layer:
geometries.AddGeometry(feature.GetGeometryRef())
if os.path.exists(output_name):
driver.DeleteDataSource(output_name)
datasource = driver.CreateDataSource(output_name)
out_layer = datasource.CreateLayer(str('states_convexhull'), spatial_reference, geom_type=ogr.wkbPolygon)
out_layer.CreateField(ogr.FieldDefn(str('id'), ogr.OFTInteger))
featureDefn = out_layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(geometries.ConvexHull())
feature.SetField(str('id'), 1)
out_layer.CreateFeature(feature)
shape.Destroy()
datasource.Destroy()
def pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='', targetSR=''):
# If you want to garuntee lon lat output, specify TargetSR otherwise, geocoords will be in image geo reference
# targetSR = osr.SpatialReference()
# targetSR.ImportFromEPSG(4326)
# Transform can be performed at the polygon level instead of pixel level
if targetSR =='':
performReprojection=False
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326)
else:
performReprojection=True
if geomTransform=='':
srcRaster = gdal.Open(inputRaster)
geomTransform = srcRaster.GetGeoTransform()
source_sr = osr.SpatialReference()
source_sr.ImportFromWkt(srcRaster.GetProjectionRef())
geom = ogr.Geometry(ogr.wkbPoint)
xOrigin = geomTransform[0]
yOrigin = geomTransform[3]
pixelWidth = geomTransform[1]
pixelHeight = geomTransform[5]
xCoord = (xPix * pixelWidth) + xOrigin
yCoord = (yPix * pixelHeight) + yOrigin
geom.AddPoint(xCoord, yCoord)
if performReprojection:
if sourceSR=='':
srcRaster = gdal.Open(inputRaster)
sourceSR = osr.SpatialReference()
sourceSR.ImportFromWkt(srcRaster.GetProjectionRef())
coord_trans = osr.CoordinateTransformation(sourceSR, targetSR)
geom.Transform(coord_trans)
return (geom.GetX(), geom.GetY())
def geoPolygonToPixelPolygonWKT(geom, inputRaster, targetSR, geomTransform, breakMultiPolygonGeo=True,
pixPrecision=2):
# Returns Pixel Coordinate List and GeoCoordinateList
polygonPixBufferList = []
polygonPixBufferWKTList = []
polygonGeoWKTList = []
if geom.GetGeometryName() == 'POLYGON':
polygonPix = ogr.Geometry(ogr.wkbPolygon)
for ring in geom:
# GetPoint returns a tuple not a Geometry
ringPix = ogr.Geometry(ogr.wkbLinearRing)
for pIdx in xrange(ring.GetPointCount()):
lon, lat, z = ring.GetPoint(pIdx)
xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)
xPix = round(xPix, pixPrecision)
yPix = round(yPix, pixPrecision)
ringPix.AddPoint(xPix, yPix)
polygonPix.AddGeometry(ringPix)
polygonPixBuffer = polygonPix.Buffer(0.0)
polygonPixBufferList.append([polygonPixBuffer, geom])
elif geom.GetGeometryName() == 'MULTIPOLYGON':
for poly in geom:
polygonPix = ogr.Geometry(ogr.wkbPolygon)
for ring in poly:
# GetPoint returns a tuple not a Geometry
ringPix = ogr.Geometry(ogr.wkbLinearRing)
for pIdx in xrange(ring.GetPointCount()):
lon, lat, z = ring.GetPoint(pIdx)
xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)
xPix = round(xPix, pixPrecision)
yPix = round(yPix, pixPrecision)
ringPix.AddPoint(xPix, yPix)
polygonPix.AddGeometry(ringPix)
polygonPixBuffer = polygonPix.Buffer(0.0)
if breakMultiPolygonGeo:
polygonPixBufferList.append([polygonPixBuffer, poly])
else:
polygonPixBufferList.append([polygonPixBuffer, geom])
for polygonTest in polygonPixBufferList:
if polygonTest[0].GetGeometryName() == 'POLYGON':
polygonPixBufferWKTList.append([polygonTest[0].ExportToWkt(), polygonTest[1].ExportToWkt()])
elif polygonTest[0].GetGeometryName() == 'MULTIPOLYGON':
for polygonTest2 in polygonTest[0]:
polygonPixBufferWKTList.append([polygonTest2.ExportToWkt(), polygonTest[1].ExportToWkt()])
return polygonPixBufferWKTList
def polygon_to_shapefile(polygons, poly_sr, fname, fields_defs=None,
poly_attrs=None):
"""
Write a set of polygons to a shapefile
Args:
polygons: a list of (lat, lng) tuples
fields_defs: The list of fields for those polygons, as a tuple
(name, ogr type) for each field. For example :
[('field1', ogr.OFTReal), ('field2', ogr.OFTInteger)]
poly_attrs: A list of dict containing the attributes for each polygon
[{'field1' : 1.0, 'field2': 42},
{'field1' : 3.0, 'field2': 60}]
"""
shp_driver = ogr.GetDriverByName("ESRI Shapefile")
out_ds = shp_driver.CreateDataSource(fname)
assert out_ds is not None, "Failed to create temporary %s" % fname
out_layer = out_ds.CreateLayer(fname, poly_sr, geom_type=ogr.wkbPolygon)
has_attrs = fields_defs is not None
if has_attrs:
attrs_name = []
for field_name, field_type in fields_defs:
out_layer.CreateField(ogr.FieldDefn(field_name, field_type))
attrs_name.append(field_name)
layer_defn = out_layer.GetLayerDefn()
for i, poly in enumerate(polygons):
ring = ogr.Geometry(ogr.wkbLinearRing)
# gdal uses the (x, y) convention => (lng, lat)
for point in poly:
ring.AddPoint(point[1], point[0])
ring.AddPoint(poly[0][1], poly[0][0]) # re-add the start to close
p = ogr.Geometry(ogr.wkbPolygon)
p.AddGeometry(ring)
out_feature = ogr.Feature(layer_defn)
out_feature.SetGeometry(p)
if has_attrs:
attrs = poly_attrs[i]
for field_name in attrs_name:
out_feature.SetField(field_name, attrs[field_name])
out_layer.CreateFeature(out_feature)
out_feature.Destroy()
out_ds.Destroy()
def filterByCoverage(vectorFile, rasterFile, covPerc):
srcVector = ogr.Open(vectorFile)
srcLayer = srcVector.GetLayer()
# merge all features in one geometry (multi polygone)
multi = ogr.Geometry(ogr.wkbMultiPolygon)
for feature in srcLayer:
geom = feature.GetGeometryRef()
multi.AddGeometry(geom)
#attributes of raster file
rasList = raster2array(rasterFile)
xsize = rasList[4][0]
ysize = abs(rasList[4][1])
pixel_area = xsize*ysize
rows = rasList[0].shape[0]
cols = rasList[0].shape[1]
x1 = rasList[2][0]
y1 = rasList[2][3]
#iterate over raster cells
for i in range(rows):
for j in range(cols):
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(x1, y1)
ring.AddPoint(x1 + xsize, y1)
ring.AddPoint(x1 + xsize, y1 - ysize)
ring.AddPoint(x1, y1 - ysize)
ring.AddPoint(x1, y1)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
intersect = multi.Intersection(poly)
if intersect.ExportToWkt() != 'GEOMETRYCOLLECTION EMPTY':
perc = (intersect.GetArea()/pixel_area)*100
if perc > covPerc:
rasList[0][i][j] = numpy.nan
x1 += xsize
x1 = rasList[2][0]
y1 -= ysize
return(rasList[0]) #return the filtered array
# numpy array to geo raster
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
poGeom = None
poDS = ogr.Open( pszDS, False )
if poDS is None:
return None
if pszSQL is not None:
poLyr = poDS.ExecuteSQL( pszSQL, None, None )
elif pszLyr is not None:
poLyr = poDS.GetLayerByName(pszLyr)
else:
poLyr = poDS.GetLayer(0)
if poLyr is None:
print("Failed to identify source layer from datasource.")
poDS.Destroy()
return None
if pszWhere is not None:
poLyr.SetAttributeFilter(pszWhere)
poFeat = poLyr.GetNextFeature()
while poFeat is not None:
poSrcGeom = poFeat.GetGeometryRef()
if poSrcGeom is not None:
eType = wkbFlatten(poSrcGeom.GetGeometryType())
if poGeom is None:
poGeom = ogr.Geometry( ogr.wkbMultiPolygon )
if eType == ogr.wkbPolygon:
poGeom.AddGeometry( poSrcGeom )
elif eType == ogr.wkbMultiPolygon:
for iGeom in range(poSrcGeom.GetGeometryCount()):
poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )
else:
print("ERROR: Geometry not of polygon type." )
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return None
poFeat = poLyr.GetNextFeature()
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return poGeom
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
poGeom = None
poDS = ogr.Open( pszDS, False )
if poDS is None:
return None
if pszSQL is not None:
poLyr = poDS.ExecuteSQL( pszSQL, None, None )
elif pszLyr is not None:
poLyr = poDS.GetLayerByName(pszLyr)
else:
poLyr = poDS.GetLayer(0)
if poLyr is None:
print("Failed to identify source layer from datasource.")
poDS.Destroy()
return None
if pszWhere is not None:
poLyr.SetAttributeFilter(pszWhere)
poFeat = poLyr.GetNextFeature()
while poFeat is not None:
poSrcGeom = poFeat.GetGeometryRef()
if poSrcGeom is not None:
eType = wkbFlatten(poSrcGeom.GetGeometryType())
if poGeom is None:
poGeom = ogr.Geometry( ogr.wkbMultiPolygon )
if eType == ogr.wkbPolygon:
poGeom.AddGeometry( poSrcGeom )
elif eType == ogr.wkbMultiPolygon:
for iGeom in range(poSrcGeom.GetGeometryCount()):
poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )
else:
print("ERROR: Geometry not of polygon type." )
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return None
poFeat = poLyr.GetNextFeature()
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return poGeom
def _create_dst_features(self, dst, trg, **kwargs):
""" Create needed OGR.Features in dst OGR.Layer
Parameters
----------
dst : OGR.Layer
destination layer
trg : OGR.Geometry
target polygon
"""
# TODO: kwargs necessary?
# claim and reset source ogr layer
layer = self.src.ds.GetLayerByName('src')
layer.ResetReading()
# if given, we apply a buffer value to the target polygon filter
trg_index = trg.GetField('index')
trg = trg.GetGeometryRef()
trg = trg.Buffer(self._buffer)
layer.SetSpatialFilter(trg)
# iterate over layer features
for ogr_src in layer:
geom = ogr_src.GetGeometryRef()
# calculate intersection, if not fully contained
if not trg.Contains(geom):
geom = trg.Intersection(geom)
# checking GeometryCollection, convert to only Polygons,
# Multipolygons
if geom.GetGeometryType() in [7]:
geocol = georef.ogr_geocol_to_numpy(geom)
geom = georef.numpy_to_ogr(geocol, 'MultiPolygon')
# only geometries containing points
if geom.IsEmpty():
continue
if geom.GetGeometryType() in [3, 6, 12]:
idx = ogr_src.GetField('index')
georef.ogr_add_geometry(dst, geom, [idx, trg_index])
def crop(seq, maxpoints=20, proximity=100, straighten=False):
x = map(float, range(0, len(seq)))
VWpts = simplify(x, seq, maxpoints)
xn, yn = map(list, zip(*VWpts))
simple = np.interp(x, xn, yn)
seq[abs(seq - simple) > proximity] = simple[abs(seq - simple) > proximity]
points = []
for xi, yi in enumerate(seq):
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(xi, yi)
points.append(point)
points = np.array(points)
while True:
poly = createPoly(xn, yn, len(seq), max(seq))
line = ogr.Geometry(ogr.wkbLineString)
for xi, yi in zip(xn, yn):
line.AddPoint(xi, yi)
dists = np.array([line.Distance(point) for point in points])
contain = np.array([point.Within(poly) for point in points])
dists[~contain] = 0
points = points[(dists > 0)]
dists = dists[(dists > 0)]
if len(dists) == 0:
break
candidate = points[np.argmax(dists)]
cp = candidate.GetPoint()
index = np.argmin(np.array(xn) < cp[0])
xn.insert(index, cp[0])
yn.insert(index, cp[1])
if straighten:
indices = [i for i in range(0, len(xn)) if (xn[i], yn[i]) in VWpts]
for i, j in enumerate(indices):
if i < (len(indices)-1):
if indices[i+1] > j+1:
dx = abs(xn[j] - xn[indices[i + 1]])
dy = abs(yn[j] - yn[indices[i + 1]])
if dx > dy:
seg_y = yn[j:indices[i + 1]+1]
print(seg_y)
for k in range(j, indices[i + 1]+1):
yn[k] = min(seg_y)
return np.interp(x, xn, yn)
def reduce(seq, maxpoints=20, straighten=False):
if min(seq) == max(seq):
return np.array(seq)
x = map(float, range(0, len(seq)))
# plt.plot(seq)
VWpts = simplify(x, seq, maxpoints)
xn, yn = map(list, zip(*VWpts))
# plt.plot(xn, yn, linewidth=2, color='r')
simple = np.interp(x, xn, yn)
seq2 = np.copy(seq)
points = []
for xi, yi in enumerate(seq2):
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(xi, yi)
points.append(point)
points = np.array(points)
while True:
poly = createPoly(xn, yn, len(seq2), max(seq2))
line = ogr.Geometry(ogr.wkbLineString)
for xi, yi in zip(xn, yn):
line.AddPoint(xi, yi)
dists = np.array([line.Distance(point) for point in points])
contain = np.array([point.Within(poly) for point in points])
dists[~contain] = 0
points = points[(dists > 0)]
dists = dists[(dists > 0)]
if len(dists) == 0:
break
candidate = points[np.argmax(dists)]
cp = candidate.GetPoint()
index = np.argmin(np.array(xn) < cp[0])
xn.insert(index, cp[0])
yn.insert(index, cp[1])
if straighten:
indices = [i for i in range(0, len(xn)) if (xn[i], yn[i]) in VWpts]
for i, j in enumerate(indices):
if i < (len(indices)-1):
if indices[i+1] > j+1:
dx = abs(xn[j] - xn[indices[i + 1]])
dy = abs(yn[j] - yn[indices[i + 1]])
if dx > dy:
seg_y = yn[j:indices[i + 1]+1]
for k in range(j, indices[i + 1]+1):
yn[k] = min(seg_y)
# plt.plot(xn, yn, linewidth=2, color='g')
# plt.show()
return np.interp(x, xn, yn).astype(int)
def point_to_pixel_geometry(points, source_epsg=None, target_epsg=None, pixel_side_length=30):
'''
Where points is a list of X,Y tuples and X and Y are coordinates in
meters, returns a series of OGR Polygons where each Polygon is the
pixel extent with a given point at its center. Assumes square pixels.
Arguments:
points Sequence of X,Y numeric pairs or OGR Point geometries
source_epsg The EPSG code of the source projection (Optional)
target_epsg The EPSG code of the target projection (Optional)
pixel_side_length The length of one side of the intended pixel
'''
polys = []
# Convert points to atomic X,Y pairs if necessary
if isinstance(points[0], ogr.Geometry):
points = [(p.GetX(), p.GetY()) for p in points]
source_ref = target_ref = None
if all((source_epsg, target_epsg)):
source_ref = osr.SpatialReference()
target_ref = osr.SpatialReference()
source_ref.ImportFromEPSG(source_epsg)
target_ref.ImportFromEPSG(target_epsg)
transform = osr.CoordinateTransformation(source_ref, target_ref)
for p in points:
r = pixel_side_length / 2 # Half the pixel width
ring = ogr.Geometry(ogr.wkbLinearRing) # Create a ring
vertices = [
(p[0] - r, p[1] + r), # Top-left
(p[0] + r, p[1] + r), # Top-right, clockwise from here...
(p[0] + r, p[1] - r),
(p[0] - r, p[1] - r),
(p[0] - r, p[1] + r) # Add top-left again to close ring
]
# Coordinate transformation
if all((source_ref, target_ref)):
vertices = [transform.TransformPoint(*v)[0:2] for v in vertices]
for vert in vertices:
ring.AddPoint(*vert)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
polys.append(poly)
return polys
