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类Feature()的实例源码
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 write(self, outfile, format="ESRI Shapefile", overwrite=True):
(outfilepath, outfilename) = os.path.split(outfile)
basename = os.path.splitext(outfilename)[0]
driver = ogr.GetDriverByName(format)
if os.path.exists(outfile):
if overwrite:
driver.DeleteDataSource(outfile)
else:
raise IOError("file already exists")
outdataset = driver.CreateDataSource(outfile)
outlayer = outdataset.CreateLayer(self.layername, geom_type=self.geomType)
outlayerdef = outlayer.GetLayerDefn()
for fieldDef in self.fieldDefs:
outlayer.CreateField(fieldDef)
for feature in self.layer:
outFeature = ogr.Feature(outlayerdef)
outFeature.SetGeometry(feature.GetGeometryRef())
for j in range(0, self.nfields):
outFeature.SetField(self.fieldnames[j], feature.GetField(j))
# add the feature to the shapefile
outlayer.CreateFeature(outFeature)
outFeature.Destroy()
srs_out = self.srs.Clone()
srs_out.MorphToESRI()
with open(os.path.join(outfilepath, basename+".prj"), "w") as prj:
prj.write(srs_out.ExportToWkt())
outdataset.Destroy()
def get_idx_as_shp(self, path, gt, wkt):
'''
Exports a Shapefile containing the locations of the extracted
endmembers. Assumes the coordinates are in decimal degrees.
'''
coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(path)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
for pair in coords:
feature = ogr.Feature(layer.GetLayerDefn())
# Create the point from the Well Known Text
point = ogr.CreateGeometryFromWkt('POINT(%f %f)' % pair)
feature.SetGeometry(point) # Set the feature geometry
layer.CreateFeature(feature) # Create the feature in the layer
feature.Destroy() # Destroy the feature to free resources
# Destroy the data source to free resources
ds.Destroy()
def buildTindex(rasterFolder, rasterExtention='.tif'):
rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
print(rasterList)
print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
memDriver = ogr.GetDriverByName('MEMORY')
gTindex = memDriver.CreateDataSource('gTindex')
srcImage = gdal.Open(rasterList[0])
spat_ref = osr.SpatialReference()
spat_ref.SetProjection(srcImage.GetProjection())
gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)
# Add an ID field
idField = ogr.FieldDefn("location", ogr.OFTString)
gTindexLayer.CreateField(idField)
# Create the feature and set values
featureDefn = gTindexLayer.GetLayerDefn()
for rasterFile in rasterList:
srcImage = gdal.Open(rasterFile)
geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)
feature = ogr.Feature(featureDefn)
feature.SetGeometry(polyToCut)
feature.SetField("location", rasterFile)
gTindexLayer.CreateFeature(feature)
feature = None
return gTindex, gTindexLayer
def ogr_add_feature(ds, src, name=None):
""" Creates OGR.Feature objects in OGR.Layer object.
.. versionadded:: 0.7.0
OGR.Features are built from numpy src points or polygons.
OGR.Features 'FID' and 'index' corresponds to source data element
Parameters
----------
ds : gdal.Dataset
object
src : :func:`numpy:numpy.array`
source data
name : string
name of wanted Layer
"""
if name is not None:
lyr = ds.GetLayerByName(name)
else:
lyr = ds.GetLayer()
defn = lyr.GetLayerDefn()
geom_name = ogr.GeometryTypeToName(lyr.GetGeomType())
fields = [defn.GetFieldDefn(i).GetName()
for i in range(defn.GetFieldCount())]
feat = ogr.Feature(defn)
for index, src_item in enumerate(src):
geom = numpy_to_ogr(src_item, geom_name)
if 'index' in fields:
feat.SetField('index', index)
feat.SetGeometry(geom)
lyr.CreateFeature(feat)
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 addfeature(self, geometry, fieldname, fieldvalue):
if fieldname not in self.fieldnames:
raise IOError("field does not exist")
featureDefn = self.layerdef
feature = ogr.Feature(featureDefn)
feature.SetGeometry(geometry)
feature.SetField(fieldname, fieldvalue)
self.layer.CreateFeature(feature)
feature.Destroy()
self.init_features()
def export_shp_nodes(self, outshp):
from osgeo import ogr
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(outshp)
layer = ds.CreateLayer('', None, ogr.wkbPoint25D)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
## If there are multiple geometries, put the "for" loop here
for id, point in self.tin_nodes.iteritems():
# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', int(id))
# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(point.wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
# Save and close everything
ds = layer = feat = geom = None
return
def export_shp_triangles(self, outshp):
from osgeo import ogr
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(outshp)
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
## If there are multiple geometries, put the "for" loop here
for id, poly in self.triangles.iteritems():
# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', int(id))
# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(poly.wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
# Save and close everything
ds = layer = feat = geom = None
def export_shp_hull(self, outshp):
from osgeo import ogr
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(outshp)
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()
## If there are multiple geometries, put the "for" loop here
for id, poly in self.hull_polygons.iteritems():
# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', int(id))
# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(poly.wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
# Save and close everything
ds = layer = feat = geom = None
def export_shp_breaklines(self, outshp):
from osgeo import ogr
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(outshp)
layer = ds.CreateLayer('', None, ogr.wkbLineString25D)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
fieldLineType = ogr.FieldDefn('LineType', ogr.OFTString)
fieldLineType.SetWidth(4)
layer.CreateField(fieldLineType)
defn = layer.GetLayerDefn()
## If there are multiple geometries, put the "for" loop here
for id, line in self.breaklines.iteritems():
# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', int(id))
feat.SetField('LineType', line['linetype'])
# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(line['geometry'].wkb)
feat.SetGeometry(geom)
layer.CreateFeature(feat)
feat = geom = None # destroy these
# Save and close everything
ds = layer = feat = geom = None
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 copy_features(src_layer, dst_layer, fix_geometry, simplify_geometry, start_index, total):
index = 0
batch_size = 200
index_batch = 0
for feat in src_layer:
if index < start_index:
index = index + 1
continue
try:
geom = shapely.wkt.loads(feat.GetGeometryRef().ExportToWkt())
except Exception as e:
print('Error({0}), skipping geometry.'.format(e))
continue
if fix_geometry and not geom.is_valid:
geom = geom.buffer(0.0)
if simplify_geometry:
geom = geom.simplify(0.004)
f = ogr.Feature(dst_layer.GetLayerDefn())
# set field values
for i in range(feat.GetFieldCount()):
fd = feat.GetFieldDefnRef(i)
f.SetField(fd.GetName(), feat.GetField(fd.GetName()))
# set geometry
f.SetGeometry(ogr.CreateGeometryFromWkt(geom.to_wkt()))
if index_batch == 0:
dst_layer.StartTransaction()
# create feature
feature = dst_layer.CreateFeature(f)
f.Destroy()
index_batch = index_batch + 1
if index_batch >= batch_size or index == total - 1:
dst_layer.CommitTransaction()
count = dst_layer.GetFeatureCount() # update number of inserted features
print('Inserted {0} of {1} features ({2:.2f}%)'.format(count, total, 100. * float(count) / total))
index_batch = 0
if index == total - 1:
break
index = index + 1
def convertLabelStringToPoly(shapeFileSrc, outGeoJSon, labelType='Airplane'):
shapeSrc = ogr.Open(shapeFileSrc)
source_layer = shapeSrc.GetLayer()
source_srs = source_layer.GetSpatialRef()
# Create the output Layer
outDriver = ogr.GetDriverByName("geojson")
if os.path.exists(outGeoJSon):
outDriver.DeleteDataSource(outGeoJSon)
outDataSource = outDriver.CreateDataSource(outGeoJSon)
outLayer = outDataSource.CreateLayer("groundTruth", source_srs, geom_type=ogr.wkbPolygon)
# Add input Layer Fields to the output Layer
inLayerDefn = source_layer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outLayer.CreateField(fieldDefn)
outLayer.CreateField(ogr.FieldDefn("Length_m", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("Width_m", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("Aspect(L/W)", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("compassDeg", ogr.OFTReal))
outLayerDefn = outLayer.GetLayerDefn()
for inFeature in source_layer:
outFeature = ogr.Feature(outLayerDefn)
for i in range(0, inLayerDefn.GetFieldCount()):
outFeature.SetField(inLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
geom = inFeature.GetGeometryRef()
if labelType == 'Airplane':
poly, Length, Width, Aspect, Direction = evaluateLineStringPlane(geom, label='Airplane')
elif labelType == 'Boat':
poly, Length, Width, Aspect, Direction = evaluateLineStringBoat(geom, label='Boat')
outFeature.SetGeometry(poly)
outFeature.SetField("Length_m", Length)
outFeature.SetField("Width_m", Width)
outFeature.SetField("Aspect(L/W)", Aspect)
outFeature.SetField("compassDeg", Direction)
outLayer.CreateFeature(outFeature)
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 lyr_proj(lyr, t_srs, preserve_fields=True):
"""Reproject an OGR layer
"""
#Need to check t_srs
s_srs = lyr.GetSpatialRef()
cT = osr.CoordinateTransformation(s_srs, t_srs)
#Do everything in memory
drv = ogr.GetDriverByName('Memory')
#Might want to save clipped, warped shp to disk?
# create the output layer
#drv = ogr.GetDriverByName('ESRI Shapefile')
#out_fn = '/tmp/temp.shp'
#if os.path.exists(out_fn):
# driver.DeleteDataSource(out_fn)
#out_ds = driver.CreateDataSource(out_fn)
out_ds = drv.CreateDataSource('out')
outlyr = out_ds.CreateLayer('out', srs=t_srs, geom_type=lyr.GetGeomType())
if preserve_fields:
# add fields
inLayerDefn = lyr.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outlyr.CreateField(fieldDefn)
# get the output layer's feature definition
outLayerDefn = outlyr.GetLayerDefn()
# loop through the input features
inFeature = lyr.GetNextFeature()
while inFeature:
# get the input geometry
geom = inFeature.GetGeometryRef()
# reproject the geometry
geom.Transform(cT)
# create a new feature
outFeature = ogr.Feature(outLayerDefn)
# set the geometry and attribute
outFeature.SetGeometry(geom)
if preserve_fields:
for i in range(0, outLayerDefn.GetFieldCount()):
outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
# add the feature to the shapefile
outlyr.CreateFeature(outFeature)
# destroy the features and get the next input feature
inFeature = lyr.GetNextFeature()
#NOTE: have to operate on ds here rather than lyr, otherwise segfault
return out_ds
#See https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#convert-vector-layer-to-array
#Should check srs, as shp could be WGS84
def geom2shp(geom, out_fn, fields=False):
"""Write out a new shapefile for input geometry
"""
from pygeotools.lib import timelib
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName(driverName)
if os.path.exists(out_fn):
drv.DeleteDataSource(out_fn)
out_ds = drv.CreateDataSource(out_fn)
out_lyrname = os.path.splitext(os.path.split(out_fn)[1])[0]
geom_srs = geom.GetSpatialReference()
geom_type = geom.GetGeometryType()
out_lyr = out_ds.CreateLayer(out_lyrname, geom_srs, geom_type)
if fields:
field_defn = ogr.FieldDefn("name", ogr.OFTString)
field_defn.SetWidth(128)
out_lyr.CreateField(field_defn)
field_defn = ogr.FieldDefn("path", ogr.OFTString)
field_defn.SetWidth(254)
out_lyr.CreateField(field_defn)
#field_defn = ogr.FieldDefn("date", ogr.OFTString)
#This allows sorting by date
field_defn = ogr.FieldDefn("date", ogr.OFTInteger)
field_defn.SetWidth(32)
out_lyr.CreateField(field_defn)
field_defn = ogr.FieldDefn("decyear", ogr.OFTReal)
field_defn.SetPrecision(8)
field_defn.SetWidth(64)
out_lyr.CreateField(field_defn)
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
out_feat.SetGeometry(geom)
if fields:
#Hack to force output extesion to tif, since out_fn is shp
out_path = os.path.splitext(out_fn)[0] + '.tif'
out_feat.SetField("name", os.path.split(out_path)[-1])
out_feat.SetField("path", out_path)
#Try to extract a date from input raster fn
out_feat_date = timelib.fn_getdatetime(out_fn)
if out_feat_date is not None:
datestamp = int(out_feat_date.strftime('%Y%m%d'))
#out_feat_date = int(out_feat_date.strftime('%Y%m%d%H%M'))
out_feat.SetField("date", datestamp)
decyear = timelib.dt2decyear(out_feat_date)
out_feat.SetField("decyear", decyear)
out_lyr.CreateFeature(out_feat)
out_ds = None
#return status?
def split_shape_into_features(shape_name, destination_directory, column, name, sep):
'''
This method will take a input shape and iterate over its features, creating
a new shape file with each one of them. It copies all the fields and the
same spatial reference from the original file. The created files are saved
in the destination directory using the number of the field given.
'''
driver = ogr.GetDriverByName(str('ESRI Shapefile'))
shape = driver.Open(shape_name, 0)
layer = shape.GetLayer()
layer_name = layer.GetName()
spatial_reference = layer.GetSpatialRef()
in_feature = layer.GetNextFeature()
shape_files = []
while in_feature:
encoding = 'utf-8'
in_feature_name = in_feature.GetField(column)
in_name = create_filename_from_string(in_feature.GetField(name).decode(encoding))
final_path = destination_directory + str(in_feature.GetField(0))
create_directory_path(final_path)
output_name = create_filename(final_path, '%s__%s.shp' % (in_feature_name, in_name))
shape_files.append(output_name)
if os.path.exists(output_name):
driver.DeleteDataSource(output_name)
datasource = driver.CreateDataSource(output_name)
outLayer = datasource.CreateLayer(layer_name, spatial_reference, geom_type=ogr.wkbPolygon)
inLayerDefn = layer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
#LOGGER.debug(fieldDefn.GetName())
outLayer.CreateField(fieldDefn)
outLayerDefn = outLayer.GetLayerDefn()
geometry = in_feature.GetGeometryRef()
out_feature = ogr.Feature(outLayerDefn)
out_feature.SetGeometry(geometry)
for i in range(0, outLayerDefn.GetFieldCount()):
out_feature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))
outLayer.CreateFeature(out_feature)
out_feature.Destroy()
in_feature.Destroy()
in_feature = layer.GetNextFeature()
shape.Destroy()
datasource.Destroy()
return shape_files
