def createGeoJSONFromRaster(geoJsonFileName, array2d, geom, proj,
layerName="BuildingID",
fieldName="BuildingID"):
memdrv = gdal.GetDriverByName('MEM')
src_ds = memdrv.Create('', array2d.shape[1], array2d.shape[0], 1)
src_ds.SetGeoTransform(geom)
src_ds.SetProjection(proj)
band = src_ds.GetRasterBand(1)
band.WriteArray(array2d)
dst_layername = "BuildingID"
drv = ogr.GetDriverByName("geojson")
dst_ds = drv.CreateDataSource(geoJsonFileName)
dst_layer = dst_ds.CreateLayer(layerName, srs=None)
fd = ogr.FieldDefn(fieldName, ogr.OFTInteger)
dst_layer.CreateField(fd)
dst_field = 1
gdal.Polygonize(band, None, dst_layer, dst_field, [], callback=None)
return
python类FieldDefn()的实例源码
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
def set_attribute(self, name, values):
""" Add/Set given Attribute with given values
Parameters
----------
name : string
Attribute Name
values : :class:`numpy:numpy.ndarray`
Values to fill in attributes
"""
lyr = self.ds.GetLayerByIndex(0)
lyr.ResetReading()
# todo: automatically check for value type
defn = lyr.GetLayerDefn()
if defn.GetFieldIndex(name) == -1:
lyr.CreateField(ogr.FieldDefn(name, ogr.OFTReal))
for i, item in enumerate(lyr):
item.SetField(name, values[i])
lyr.SetFeature(item)
def polygonize(self,rasterTemp,outShp):
sourceRaster = gdal.Open(rasterTemp)
band = sourceRaster.GetRasterBand(1)
driver = ogr.GetDriverByName("ESRI Shapefile")
# If shapefile already exist, delete it
if os.path.exists(outShp):
driver.DeleteDataSource(outShp)
outDatasource = driver.CreateDataSource(outShp)
# get proj from raster
srs = osr.SpatialReference()
srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
# create layer with proj
outLayer = outDatasource.CreateLayer(outShp,srs)
# Add class column (1,2...) to shapefile
newField = ogr.FieldDefn('Class', ogr.OFTInteger)
outLayer.CreateField(newField)
gdal.Polygonize(band, None,outLayer, 0,[],callback=None)
outDatasource.Destroy()
sourceRaster=None
band=None
ioShpFile = ogr.Open(outShp, update = 1)
lyr = ioShpFile.GetLayerByIndex(0)
lyr.ResetReading()
for i in lyr:
lyr.SetFeature(i)
# if area is less than inMinSize or if it isn't forest, remove polygon
if i.GetField('Class')!=1:
lyr.DeleteFeature(i.GetFID())
ioShpFile.Destroy()
return outShp
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_create_layer(ds, name, srs=None, geom_type=None, fields=None):
"""Creates OGR.Layer objects in gdal.Dataset object.
.. versionadded:: 0.7.0
Creates one OGR.Layer with given name in given gdal.Dataset object
using given OGR.GeometryType and FieldDefinitions
Parameters
----------
ds : gdal.Dataset
object
name : string
OGRLayer name
srs : OSR.SpatialReference
object
geom_type : OGR GeometryType
(eg. ogr.wkbPolygon)
fields : list of 2 element tuples
(strings, OGR.DataType) field name, field type
Returns
-------
out : OGR.Layer
object
"""
if geom_type is None:
raise TypeError("geometry_type needed")
lyr = ds.CreateLayer(name, srs=srs, geom_type=geom_type)
if fields is not None:
for fname, fvalue in fields:
lyr.CreateField(ogr.FieldDefn(fname, fvalue))
return lyr
def saveToShape(self,array,srs,outShapeFile):
# Parse a delimited text file of volcano data and create a shapefile
# use a dictionary reader so we can access by field name
# set up the shapefile driver
outDriver = ogr.GetDriverByName( 'ESRI Shapefile' )
# create the data source
if os.path.exists(outShapeFile):
outDriver.DeleteDataSource(outShapeFile)
# Remove output shapefile if it already exists
ds = outDriver.CreateDataSource(outShapeFile) #options = ['SPATIALITE=YES'])
# create the spatial reference, WGS84
lyrout = ds.CreateLayer('randomSubset',srs)
fields = [array[1].GetFieldDefnRef(i).GetName() for i in range(array[1].GetFieldCount())]
for f in fields:
field_name = ogr.FieldDefn(f, ogr.OFTString)
field_name.SetWidth(24)
lyrout.CreateField(field_name)
for k in array:
lyrout.CreateFeature(k)
# Save and close the data source
ds = None
def addfield(self, name, type=ogr.OFTString, width=10):
fieldDefn = ogr.FieldDefn(name, type)
if type == ogr.OFTString:
fieldDefn.SetWidth(width)
self.layer.CreateField(fieldDefn)
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 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 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 _create_layer(self, type, srid, attributes, title):
"""
Creates an empty spatialite layer
:param type: 'Point', 'LineString', 'Polygon', etc.
:param srid: CRS ID of the layer
:param attributes: list of (attribute_name, attribute_type, attribute_typename)
:param title: title of the layer
"""
driver = ogr.GetDriverByName('GPKG')
ds = driver.CreateDataSource(self.output_local_file)
layer = ds.CreateLayer("meta", geom_type = ogr.wkbNone)
layer.CreateField(ogr.FieldDefn('key', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('value', ogr.OFTString))
if srid:
wkbType = { 'point': ogr.wkbPoint,
'multipoint': ogr.wkbMultiPoint,
'linestring': ogr.wkbLineString,
'multilinestring': ogr.wkbMultiLineString,
'polygon': ogr.wkbPolygon,
'multipolygon': ogr.wkbMultiPolygon
}[type]
srs = osr.SpatialReference()
srs.ImportFromEPSGA(int(srid))
else:
wkbType = ogr.wkbNone
srs = None
layer = ds.CreateLayer("data", srs, wkbType, ['FID=id'])
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64))
layer.CreateField(ogr.FieldDefn('fid', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('_xml_', ogr.OFTString))
att_type_map = {QVariant.String : ogr.OFTString,
QVariant.Int : ogr.OFTInteger,
QVariant.Double: ogr.OFTReal,
QVariant.DateTime: ogr.OFTDateTime}
for aname, atype in attributes:
layer.CreateField(ogr.FieldDefn(aname, att_type_map[atype]))
# update fields
layer.ResetReading()
qgs_layer = QgsVectorLayer("{}|layername=data".format(self.output_local_file), title, "ogr")
return qgs_layer