def exporttogeojson(geojsonfilename, buildinglist):
#
# geojsonname should end with .geojson
# building list should be list of dictionaries
# list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'polyPix': poly,
# 'polyGeo': poly}
# image_id is a string,
# BuildingId is an integer,
# poly is a ogr.Geometry Polygon
#
# returns geojsonfilename
driver = ogr.GetDriverByName('geojson')
if os.path.exists(geojsonfilename):
driver.DeleteDataSource(geojsonfilename)
datasource = driver.CreateDataSource(geojsonfilename)
layer = datasource.CreateLayer('buildings', geom_type=ogr.wkbPolygon)
field_name = ogr.FieldDefn("ImageId", ogr.OFTString)
field_name.SetWidth(75)
layer.CreateField(field_name)
layer.CreateField(ogr.FieldDefn("BuildingId", ogr.OFTInteger))
# loop through buildings
for building in buildinglist:
# create feature
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField("ImageId", building['ImageId'])
feature.SetField("BuildingId", building['BuildingId'])
feature.SetGeometry(building['polyPix'])
# Create the feature in the layer (geojson)
layer.CreateFeature(feature)
# Destroy the feature to free resources
feature.Destroy()
datasource.Destroy()
return geojsonfilename
python类Geometry()的实例源码
def mergePolyList(geojsonfilename):
multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
datasource = ogr.Open(geojsonfilename, 0)
layer = datasource.GetLayer()
print(layer.GetFeatureCount())
for idx, feature in enumerate(layer):
poly = feature.GetGeometryRef()
if poly:
multipolygon.AddGeometry(feature.GetGeometryRef().Clone())
return multipolygon
def get_envelope(poly):
env = poly.GetEnvelope()
# Get Envelope returns a tuple (minX, maxX, minY, maxY)
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(env[0], env[2])
ring.AddPoint(env[0], env[3])
ring.AddPoint(env[1], env[3])
ring.AddPoint(env[1], env[2])
ring.AddPoint(env[0], env[2])
# Create polygon
poly1 = ogr.Geometry(ogr.wkbPolygon)
poly1.AddGeometry(ring)
return poly1
def ogr_add_geometry(layer, geom, attrs):
""" Copies single OGR.Geometry object to an OGR.Layer object.
.. versionadded:: 0.7.0
Given OGR.Geometry is copied to new OGR.Feature and
written to given OGR.Layer by given index. Attributes are attached.
Parameters
----------
layer : OGR.Layer
object
geom : OGR.Geometry
object
attrs : list
attributes referring to layer fields
"""
defn = layer.GetLayerDefn()
feat = ogr.Feature(defn)
for i, item in enumerate(attrs):
feat.SetField(i, item)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
def ogr_to_numpy(ogrobj):
"""Backconvert a gdal/ogr geometry to a numpy vertex array.
.. versionadded:: 0.7.0
Using JSON as a vehicle to efficiently deal with numpy arrays.
Parameters
----------
ogrobj : ogr.Geometry
object
Returns
-------
out : :class:`numpy:numpy.ndarray`
a nested ndarray of vertices of shape (num vertices, 2)
"""
jsonobj = eval(ogrobj.ExportToJson())
return np.squeeze(jsonobj['coordinates'])
def get_centroid(polyg):
"""Return centroid of a polygon
Parameters
----------
polyg : :class:`numpy:numpy.ndarray`
of shape (num vertices, 2) or ogr.Geometry object
Returns
-------
out : x and y coordinate of the centroid
"""
if not type(polyg) == ogr.Geometry:
polyg = numpy_to_ogr(polyg, 'Polygon')
return polyg.Centroid().GetPoint()[0:2]
def get_features_as_geojson(layer, bbox=None, union=False):
""" Get features in this layer and return as GeoJSON """
if bbox is not None:
layer.SetSpatialFilterRect(bbox[0], bbox[3], bbox[2], bbox[1])
poly = ogr.Geometry(ogr.wkbPolygon)
if union:
for feature in layer:
geom = feature.GetGeometryRef()
# required for ogr2
if hasattr(geom, 'GetLinearGeometry'):
geom = geom.GetLinearGeometry()
poly = poly.Union(geom)
if bbox is not None:
wkt = "POLYGON ((%s %s, %s %s, %s %s, %s %s, %s %s))" % \
(bbox[0], bbox[1], bbox[2], bbox[1], bbox[2], bbox[3], bbox[0], bbox[3], bbox[0], bbox[1])
bbox_wkt = ogr.CreateGeometryFromWkt(wkt)
poly = poly.Intersection(bbox_wkt)
return json.loads(poly.ExportToJson())
def addMetricDistanceToEdge(x1,y1,x2,y2, epsgCode):
# we assume initial epsg is wsg84 (merctor projection)
if epsgCode != 4326:
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(4326)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(x1,y1)
line.AddPoint(x2,y2)
if epsgCode != 4326:
line.Transform(coordTrans)
length = line.Length()
return length
def addMetricDistanceToEdges(graph, epsgCode):
# we assume initial epsg is wsg84 (merctor projection)
metricDistance = {}
if epsgCode != 4326:
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(4326)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
for edge in graph.edges():
node1, node2 = edge
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
if epsgCode != 4326:
line.Transform(coordTrans)
length = line.Length()
metricDistance[edge] = length
nx.set_edge_attributes(graph, 'length', metricDistance)
return graph
def addMetricDistanceToEdges(graph, epsgCode):
# we assume initial epsg is wsg84 (merctor projection)
metricDistance = {}
if epsgCode != 4326:
sourceProjection = osr.SpatialReference()
sourceProjection.ImportFromEPSG(4326)
destinationProjection = osr.SpatialReference()
destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
for edge in graph.edges():
node1, node2 = edge
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
if epsgCode != 4326:
line.Transform(coordTrans)
length = line.Length()
metricDistance[edge] = length
nx.set_edge_attributes(graph, 'length', metricDistance)
return graph
def get_bbox(shapefile):
""" Gets the bounding box of a shapefile in EPSG 4326.
If shapefile is not in WGS84, bounds are reprojected.
"""
driver = ogr.GetDriverByName('ESRI Shapefile')
# 0 means read, 1 means write
data_source = driver.Open(shapefile, 0)
# Check to see if shapefile is found.
if data_source is None:
failure('Could not open %s' % (shapefile))
else:
layer = data_source.GetLayer()
shape_bbox = layer.GetExtent()
spatialRef = layer.GetSpatialRef()
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
# this check for non-WGS84 projections gets some false positives, but that's okay
if target.ExportToProj4() != spatialRef.ExportToProj4():
transform = osr.CoordinateTransformation(spatialRef, target)
point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(shape_bbox[0], shape_bbox[2])
point2 = ogr.Geometry(ogr.wkbPoint)
point2.AddPoint(shape_bbox[1], shape_bbox[3])
point1.Transform(transform)
point2.Transform(transform)
shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]
return shape_bbox
def readwktcsv(csv_path,removeNoBuildings=True, groundTruthFile=True):
#
# csv Format Expected = ['ImageId', 'BuildingId', 'PolygonWKT_Pix', 'PolygonWKT_Geo']
# returns list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'poly': poly}
# image_id is a string,
# BuildingId is an integer,
# poly is a ogr.Geometry Polygon
buildinglist = []
with open(csv_path, 'rb') as csvfile:
building_reader = csv.reader(csvfile, delimiter=',', quotechar='"')
next(building_reader, None) # skip the headers
for row in building_reader:
if removeNoBuildings:
if int(row[1]) != -1:
polyPix = ogr.CreateGeometryFromWkt(row[2])
if groundTruthFile:
polyGeo = ogr.CreateGeometryFromWkt(row[3])
else:
polyGeo = []
buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
'polyGeo': polyGeo,
})
else:
polyPix = ogr.CreateGeometryFromWkt(row[2])
if groundTruthFile:
polyGeo = ogr.CreateGeometryFromWkt(row[3])
else:
polyGeo = []
buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
'polyGeo': polyGeo,
})
return buildinglist
def latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
# type: (object, object, object, object, object) -> object
sourcesr = osr.SpatialReference()
sourcesr.ImportFromEPSG(4326)
geom = ogr.Geometry(ogr.wkbPoint)
geom.AddPoint(lon, lat)
if targetsr == '':
src_raster = gdal.Open(input_raster)
targetsr = osr.SpatialReference()
targetsr.ImportFromWkt(src_raster.GetProjectionRef())
coord_trans = osr.CoordinateTransformation(sourcesr, targetsr)
if geom_transform == '':
src_raster = gdal.Open(input_raster)
transform = src_raster.GetGeoTransform()
else:
transform = geom_transform
x_origin = transform[0]
# print(x_origin)
y_origin = transform[3]
# print(y_origin)
pixel_width = transform[1]
# print(pixel_width)
pixel_height = transform[5]
# print(pixel_height)
geom.Transform(coord_trans)
# print(geom.GetPoint())
x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width
y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height
return (x_pix, y_pix)
def returnBoundBox(xOff, yOff, pixDim, inputRaster, targetSR='', pixelSpace=False):
# Returns Polygon for a specific square defined by a center Pixel and
# number of pixels in each dimension.
# Leave targetSR as empty string '' or specify it as a osr.SpatialReference()
# targetSR = osr.SpatialReference()
# targetSR.ImportFromEPSG(4326)
if targetSR == '':
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326)
xCord = [xOff - pixDim / 2, xOff - pixDim / 2, xOff + pixDim / 2,
xOff + pixDim / 2, xOff - pixDim / 2]
yCord = [yOff - pixDim / 2, yOff + pixDim / 2, yOff + pixDim / 2,
yOff - pixDim / 2, yOff - pixDim / 2]
ring = ogr.Geometry(ogr.wkbLinearRing)
for idx in xrange(len(xCord)):
if pixelSpace == False:
geom = pixelToGeoCoord(xCord[idx], yCord[idx], inputRaster)
ring.AddPoint(geom[0], geom[1], 0)
else:
ring.AddPoint(xCord[idx], yCord[idx], 0)
poly = ogr.Geometry(ogr.wkbPolygon)
if pixelSpace == False:
poly.AssignSpatialReference(targetSR)
poly.AddGeometry(ring)
return poly
def search_rtree(test_building, index):
# input test poly ogr.Geometry and rtree index
if test_building.GetGeometryName() == 'POLYGON' or \
test_building.GetGeometryName() == 'MULTIPOLYGON':
fidlist = index.intersection(test_building.GetEnvelope())
else:
fidlist = []
return fidlist
def getRasterExtent(srcImage):
geoTrans = srcImage.GetGeoTransform()
ulX = geoTrans[0]
ulY = geoTrans[3]
xDist = geoTrans[1]
yDist = geoTrans[5]
rtnX = geoTrans[2]
rtnY = geoTrans[4]
cols = srcImage.RasterXSize
rows = srcImage.RasterYSize
lrX = ulX + xDist * cols
lrY = ulY + yDist * rows
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(lrX, lrY)
ring.AddPoint(lrX, ulY)
ring.AddPoint(ulX, ulY)
ring.AddPoint(ulX, lrY)
ring.AddPoint(lrX, lrY)
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
return geoTrans, poly, ulX, ulY, lrX, lrY
def createPolygonFromCorners(lrX,lrY,ulX, ulY):
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(lrX, lrY)
ring.AddPoint(lrX, ulY)
ring.AddPoint(ulX, ulY)
ring.AddPoint(ulX, lrY)
ring.AddPoint(lrX, lrY)
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
return poly
def evaluateLineStringPlane(geom, label='Airplane'):
ring = ogr.Geometry(ogr.wkbLinearRing)
for i in range(0, geom.GetPointCount()):
# GetPoint returns a tuple not a Geometry
pt = geom.GetPoint(i)
ring.AddPoint(pt[0], pt[1])
pt = geom.GetPoint(0)
ring.AddPoint(pt[0], pt[1])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
transform_WGS84_To_UTM, transform_UTM_To_WGS84, utm_cs = gT.createUTMTransform(geom)
geom.Transform(transform_WGS84_To_UTM)
pt0 = geom.GetPoint(0) # Tail
pt1 = geom.GetPoint(1) # Wing
pt2 = geom.GetPoint(2) # Nose
pt3 = geom.GetPoint(3) # Wing
Length = math.sqrt((pt2[0]-pt0[0])**2 + (pt2[1]-pt0[1])**2)
Width = math.sqrt((pt3[0] - pt1[0])**2 + (pt3[1] - pt1[1])**2)
Aspect = Length/Width
Direction = (math.atan2(pt2[0]-pt0[0], pt2[1]-pt0[1])*180/math.pi) % 360
geom.Transform(transform_UTM_To_WGS84)
return [poly, Length, Width, Aspect, Direction]
def __init__(self, lat, lng):
""" Coordinates are in degrees """
self.point = ogr.Geometry(ogr.wkbPoint)
self.point.AddPoint(lng, lat)
def transPoint(point, inproj, outproj):
sr_srs = osr.SpatialReference()
sr_srs.ImportFromEPSG(inproj)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(outproj)
coordTransform = osr.CoordinateTransformation(sr_srs, dst_srs)
gps_point = ogr.Geometry(ogr.wkbPoint)
gps_point.AddPoint(point[0], point[1])
gps_point.Transform(coordTransform)
x = float(gps_point.ExportToWkt().split()[1].split('(')[1])
y = float(gps_point.ExportToWkt().split()[2])
return [x,y]
def get_data_by_geom(self, geom=None):
""" Returns DataSource geometries filtered by given OGR geometry
Parameters
----------
geom : OGR.Geometry object
"""
lyr = self.ds.GetLayer()
lyr.ResetReading()
lyr.SetAttributeFilter(None)
lyr.SetSpatialFilter(geom)
return self._get_data()
def _get_intersection(self, trg=None, idx=None, buf=0.):
"""Just a toy function if you want to inspect the intersection
points/polygons of an arbitrary target or an target by index.
"""
# TODO: kwargs necessary?
# check wether idx is given
if idx is not None:
if self.trg:
try:
lyr = self.trg.ds.GetLayerByName('trg')
feat = lyr.GetFeature(idx)
trg = feat.GetGeometryRef()
except Exception:
raise TypeError("No target polygon found at index {0}".
format(idx))
else:
raise TypeError('No target polygons found in object!')
# check for trg
if trg is None:
raise TypeError('Either *trg* or *idx* keywords must be given!')
# check for geometry
if not type(trg) == ogr.Geometry:
trg = georef.numpy_to_ogr(trg, 'Polygon')
# apply Buffer value
trg = trg.Buffer(buf)
if idx is None:
intersecs = self.src.get_data_by_geom(trg)
else:
intersecs = self.dst.get_data_by_att('trg_index', idx)
return intersecs
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)
feat_cnt = layer.GetFeatureCount()
if feat_cnt:
[georef.ogr_add_geometry(dst, ogr_src.GetGeometryRef(),
[ogr_src.GetField('index'), trg_index])
for ogr_src in layer]
else:
layer.SetSpatialFilter(None)
src_pts = np.array([ogr_src.GetGeometryRef().GetPoint_2D(0)
for ogr_src in layer])
centroid = georef.get_centroid(trg)
tree = cKDTree(src_pts)
distnext, ixnext = tree.query([centroid[0], centroid[1]], k=1)
feat = layer.GetFeature(ixnext)
georef.ogr_add_geometry(dst, feat.GetGeometryRef(),
[feat.GetField('index'), trg_index])
def get_shape_points(geom):
"""
Extract coordinate points from given ogr geometry as generator object
If geometries are nested, function recurses.
.. versionadded:: 0.6.0
Parameters
----------
geom : ogr.Geometry
Returns
-------
result : generator object
expands to Nx2 dimensional nested point arrays
"""
type = geom.GetGeometryType()
if type:
# 1D Geometries, LINESTRINGS
if type == 2:
result = np.array(geom.GetPoints())
yield result
# RINGS, POLYGONS, MULTIPOLYGONS, MULTILINESTRINGS
elif type > 2:
# iterate over geometries and recurse
for item in geom:
for result in get_shape_points(item):
yield result
else:
print("Unknown Geometry")
def numpy_to_ogr(vert, geom_name):
"""Convert a vertex array to gdal/ogr geometry.
.. versionadded:: 0.7.0
Using JSON as a vehicle to efficiently deal with numpy arrays.
Parameters
----------
vert : array_like
a numpy array of vertices of shape (num vertices, 2)
geom_name : string
Name of Geometry
Returns
-------
out : ogr.Geometry
object of type geom_name
"""
if geom_name in ['Polygon', 'MultiPolygon']:
json_str = "{{'type':{0!r},'coordinates':[{1!r}]}}".\
format(geom_name, vert.tolist())
else:
json_str = "{{'type':{0!r},'coordinates':{1!r}}}".\
format(geom_name, vert.tolist())
return ogr.CreateGeometryFromJson(json_str)
def ogr_geocol_to_numpy(ogrobj):
"""Backconvert a gdal/ogr geometry Collection to a numpy vertex array.
.. versionadded:: 0.7.0
This extracts only Polygon geometries!
Using JSON as a vehicle to efficiently deal with numpy arrays.
Parameters
----------
ogrobj : ogr.Geometry
Collection object
Returns
-------
out : :class:`numpy:numpy.ndarray`
a nested ndarray of vertices of shape (num vertices, 2)
"""
jsonobj = eval(ogrobj.ExportToJson())
mpol = []
for item in jsonobj['geometries']:
if item['type'] == 'Polygon':
mpol.append(item['coordinates'])
return np.squeeze(mpol)
def testAddingFeatureUsingOgr(self):
dataSource =self._getOgrTestLayer()
layer = dataSource.GetLayer()
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(123, 456)
featureDefn = layer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(point)
feature.SetField("fid", 5)
feature.SetField("n", 5)
layer.CreateFeature(feature)
dataSource.Destroy()
def bbox(coordinates, crs, outname=None, format='ESRI Shapefile', overwrite=True):
"""
create a bounding box vector object or shapefile from coordinates and coordinate reference system
coordinates must be provided in a dictionary containing numerical variables with names 'xmin', 'xmax', 'ymin' and 'ymax'
the coordinate reference system can be in either WKT, EPSG or PROJ4 format
"""
srs = crsConvert(crs, 'osr')
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(coordinates['xmin'], coordinates['ymin'])
ring.AddPoint(coordinates['xmin'], coordinates['ymax'])
ring.AddPoint(coordinates['xmax'], coordinates['ymax'])
ring.AddPoint(coordinates['xmax'], coordinates['ymin'])
ring.CloseRings()
geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(ring)
geom.FlattenTo2D()
bbox = Vector(driver='Memory')
bbox.addlayer('bbox', srs, ogr.wkbPolygon)
bbox.addfield('id', width=4)
bbox.addfeature(geom, 'id', 1)
geom.Destroy()
if outname is None:
return bbox
else:
bbox.write(outname, format, overwrite)
def createPoly(xn, yn, xmax, ymax):
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(0, 0)
for item in zip(xn, yn):
item = map(int, item)
if item != [0, 0] and item != [xmax, ymax]:
ring.AddPoint(item[0], item[1])
ring.AddPoint(xmax, ymax)
ring.AddPoint(xmax, 0)
ring.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
return poly
def createPoly(xn, yn, xmax, ymax):
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(0, 0)
for item in zip(xn, yn):
item = map(int, item)
if item != [0, 0] and item != [xmax, ymax]:
ring.AddPoint(item[0], item[1])
ring.AddPoint(xmax, ymax)
ring.AddPoint(xmax, 0)
ring.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
return poly
