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()
python类GetDriverByName()的实例源码
def rasterize_shapefile_like(shpfile, model_raster_fname, dtype, options,
nodata_val=0,):
"""
Given a shapefile, rasterizes it so it has
the exact same extent as the given model_raster
`dtype` is a gdal type like gdal.GDT_Byte
`options` should be a list that will be passed to GDALRasterizeLayers
papszOptions, like ["ATTRIBUTE=vegetation","ALL_TOUCHED=TRUE"]
"""
model_dataset = gdal.Open(model_raster_fname)
shape_dataset = ogr.Open(shpfile)
shape_layer = shape_dataset.GetLayer()
mem_drv = gdal.GetDriverByName('MEM')
mem_raster = mem_drv.Create(
'',
model_dataset.RasterXSize,
model_dataset.RasterYSize,
1,
dtype
)
mem_raster.SetProjection(model_dataset.GetProjection())
mem_raster.SetGeoTransform(model_dataset.GetGeoTransform())
mem_band = mem_raster.GetRasterBand(1)
mem_band.Fill(nodata_val)
mem_band.SetNoDataValue(nodata_val)
# http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
# http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
err = gdal.RasterizeLayer(
mem_raster,
[1],
shape_layer,
None,
None,
[1],
options
)
assert err == gdal.CE_None
return mem_raster.ReadAsArray()
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 load_comboBox(self, layers):
"""Load the fields into combobox when layers are changed"""
selectedLayerIndex = self.dlg.comboBox.currentIndex()
if selectedLayerIndex < 0 or selectedLayerIndex > len(layers):
return
try:
selectedLayer = layers[selectedLayerIndex]
except:
return
fieldnames = [field.name() for field in selectedLayer.pendingFields()]
self.clear_fields()
self.dlg.comboBox_C.addItems(fieldnames)
(path, layer_id) = selectedLayer.dataProvider().dataSourceUri().split('|')
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(path, 0)
inLayer = inDataSource.GetLayer()
global type
type = inLayer.GetLayerDefn().GetGeomType()
if type == 3: # is a polygon
self.dlg.checkBox_queen.setChecked(True)
self.dlg.lineEditThreshold.setEnabled(False)
self.dlg.checkBox_optimizeDistance.setChecked(False)
self.dlg.checkBox_optimizeDistance.setEnabled(False)
self.dlg.lineEdit_minT.setEnabled(False)
self.dlg.lineEdit_maxT.setEnabled(False)
self.dlg.lineEdit_dist.setEnabled(False)
else:
self.dlg.checkBox_queen.setChecked(False)
self.dlg.lineEditThreshold.setEnabled(True)
self.dlg.checkBox_optimizeDistance.setEnabled(True)
self.dlg.lineEdit_minT.setEnabled(True)
self.dlg.lineEdit_dist.setEnabled(True)
self.dlg.lineEdit_maxT.setEnabled(True)
thresh = pysal.min_threshold_dist_from_shapefile(path)
self.dlg.lineEditThreshold.setText(str(int(thresh)))
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
def gdal_create_dataset(drv, name, cols=0, rows=0, bands=0,
gdal_type=gdal.GDT_Unknown, remove=False):
"""Creates GDAL.DataSet object.
.. versionadded:: 0.7.0
.. versionchanged:: 0.11.0
- changed parameters to keyword args
- added 'bands' as parameter
Parameters
----------
drv : string
GDAL driver string
name : string
path to filename
cols : int
# of columns
rows : int
# of rows
bands : int
# of raster bands
gdal_type : raster data type
eg. gdal.GDT_Float32
remove : bool
if True, existing gdal.Dataset will be
removed before creation
Returns
-------
out : gdal.Dataset
object
"""
driver = gdal.GetDriverByName(drv)
metadata = driver.GetMetadata()
if not metadata.get('DCAP_CREATE', False):
raise IOError("Driver %s doesn't support Create() method.".format(drv))
if remove:
if os.path.exists(name):
driver.Delete(name)
ds = driver.Create(name, cols, rows, bands, gdal_type)
return ds
def crea_gml(dxf_origen_file, gml_salida_file, src):
""" Transforma la información de la geometría de un archivo DXF al estándar de Catastro en formato GML.
:dxf_origen_file: Dirección del archivo en formato DXF con la geometría de origen
:gml_salida_file: Dirección del archivo en formato GML a sobreescribir con el resultado
:src: Sistema de Refencia de Coordendas del DXF. Según cógigos EPSG
"""
# Accede mediante GDAL al archivo DXF
driver = ogr.GetDriverByName('DXF')
data_source = driver.Open(dxf_origen_file, 0)
layer = data_source.GetLayer()
if src not in SRC_DICT: # Comprueba que el SRC es correcto
print('ERROR: El código SRC ({}) indicado es incorrecto.'.format(src))
print('Los SRC permitidos son 25828, 25829, 25830 o 25831')
sys.exit()
print('Archivo GML de salida: {}'.format(gml_salida_file))
print('Código EPSG seleccionado: {}\n'.format(src))
with open(gml_salida_file, 'w') as filegml:
filegml.writelines(PLANTILLA_1)
print('El archivo {} contiene {} geometría.'.format(dxf_origen_file, layer.GetFeatureCount()))
for feature in layer:
geom = feature.GetGeometryRef()
area = geom.Area()
print('El área del polígono es {:.4f} m2.'.format(area))
filegml.writelines(str(area)) # Añade el área al GML
perimetro = geom.GetGeometryRef(0)
print('\nTotal de vértices del polígono: {}'.format(perimetro.GetPointCount()))
print('Listado de coordenadas de los vértices:\nid,x,y')
filegml.writelines(PLANTILLA_2.format(src=src))
for i in range(0, perimetro.GetPointCount()):
pt = perimetro.GetPoint(i)
coordlist = [str(pt[0]), ' ', str(pt[1]), '\n']
filegml.writelines(coordlist) # Añade listado de coordenadas X e Y al GML
print('{},{:.4f},{:.4f}'.format(i, pt[0], pt[1]))
filegml.writelines(PLANTILLA_3) # Añade XML
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
