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类OFTInteger()的实例源码
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 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 _check_src(self, src):
""" Basic check of source elements (sequence of points or polygons).
- array cast of source elements
- create ogr_src datasource/layer holding src points/polygons
- transforming source grid points/polygons to ogr.geometries
on ogr.layer
"""
ogr_src = io.gdal_create_dataset('Memory', 'out',
gdal_type=gdal.OF_VECTOR)
try:
# is it ESRI Shapefile?
ds_in, tmp_lyr = io.open_shape(src, driver=ogr.
GetDriverByName('ESRI Shapefile'))
ogr_src_lyr = ogr_src.CopyLayer(tmp_lyr, self._name)
if self._srs is None:
self._srs = ogr_src_lyr.GetSpatialRef()
except IOError:
# no ESRI shape file
raise
# all failed? then it should be sequence or numpy array
except RuntimeError:
src = np.array(src)
# create memory datasource, layer and create features
if src.ndim == 2:
geom_type = ogr.wkbPoint
# no Polygons, just Points
else:
geom_type = ogr.wkbPolygon
fields = [('index', ogr.OFTInteger)]
georef.ogr_create_layer(ogr_src, self._name, srs=self._srs,
geom_type=geom_type, fields=fields)
georef.ogr_add_feature(ogr_src, src, name=self._name)
return ogr_src
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 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
utils.py 文件源码
项目:Python-Geospatial-Development-Third-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 19
收藏 0
点赞 0
评论 0
def get_ogr_feature_attribute(attr, feature):
attr_name = attr.name
if not feature.IsFieldSet(attr_name):
return (True, None)
if attr.type == ogr.OFTInteger:
value = str(feature.GetFieldAsInteger(attr_name))
elif attr.type == ogr.OFTIntegerList:
value = repr(feature.GetFieldAsIntegerList(attr_name))
elif attr.type == ogr.OFTReal:
value = feature.GetFieldAsDouble(attr_name)
value = "%*.*f" % (attr.width, attr.precision, value)
elif attr.type == ogr.OFTRealList:
values = feature.GetFieldAsDoubleList(attr_name)
str_values = []
for value in values:
str_values.append("%*.*f" % (attr.width,
attr.precision, value))
value = repr(str_values)
elif attr.type == ogr.OFTString:
value = feature.GetFieldAsString(attr_name)
elif attr.type == ogr.OFTStringList:
value = repr(feature.GetFieldAsStringList(attr_name))
elif attr.type == ogr.OFTDate:
parts = feature.GetFieldAsDateTime(attr_name)
year,month,day,hour,minute,second,tzone = parts
value = "%d,%d,%d,%d" % (year,month,day,tzone)
elif attr.type == ogr.OFTTime:
parts = feature.GetFieldAsDateTime(attr_name)
year,month,day,hour,minute,second,tzone = parts
value = "%d,%d,%d,%d" % (hour,minute,second,tzone)
elif attr.type == ogr.OFTDateTime:
parts = feature.GetFieldAsDateTime(attr_name)
year,month,day,hour,minute,second,tzone = parts
value = "%d,%d,%d,%d,%d,%d,%d,%d" % (year,month,day,
hour,minute,
second,tzone)
else:
return (False, "Unsupported attribute type: " +
str(attr.type))
return (True, value)
#############################################################################
utils.py 文件源码
项目:Python-Geospatial-Development-Third-Edition
作者: PacktPublishing
项目源码
文件源码
阅读 19
收藏 0
点赞 0
评论 0
def set_ogr_feature_attribute(attr, value, feature):
attr_name = attr.name
if value == None:
feature.UnsetField(attr_name)
return
if attr.type == ogr.OFTInteger:
feature.SetField(attr_name, int(value))
elif attr.type == ogr.OFTIntegerList:
integers = eval(value)
feature.SetFieldIntegerList(attr_name, integers)
elif attr.type == ogr.OFTReal:
feature.SetField(attr_name, float(value))
elif attr.type == ogr.OFTRealList:
floats = []
for s in eval(value):
floats.append(eval(s))
feature.SetFieldDoubleList(attr_name, floats)
elif attr.type == ogr.OFTString:
feature.SetField(attr_name, value)
elif attr.type == ogr.OFTStringList:
strings = []
for s in eval(value):
strings.append(s.encode(encoding))
feature.SetFieldStringList(attr_name, strings)
elif attr.type == ogr.OFTDate:
parts = value.split(",")
year = int(parts[0])
month = int(parts[1])
day = int(parts[2])
tzone = int(parts[3])
feature.SetField(attr_name, year, month, day,
0, 0, 0, tzone)
elif attr.type == ogr.OFTTime:
parts = value.split(",")
hour = int(parts[0])
minute = int(parts[1])
second = int(parts[2])
tzone = int(parts[3])
feature.SetField(attr_name, 0, 0, 0,
hour, minute, second, tzone)
elif attr.type == ogr.OFTDateTime:
parts = value.split(",")
year = int(parts[0])
month = int(parts[1])
day = int(parts[2])
hour = int(parts[3])
minute = int(parts[4])
second = int(parts[5])
tzone = int(parts[6])
feature.SetField(attr_mame, year, month, day,
hour, minute, second, tzone)