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类wkbPolygon()的实例源码
def SetZ (poGeom, dfZ ):
if poGeom is None:
return
eGType = wkbFlatten(poGeom.GetGeometryType())
if eGType == ogr.wkbPoint:
poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ)
elif eGType == ogr.wkbLineString or \
eGType == ogr.wkbLinearRing:
for i in range(poGeom.GetPointCount()):
poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ)
elif eGType == ogr.wkbPolygon or \
eGType == ogr.wkbMultiPoint or \
eGType == ogr.wkbMultiLineString or \
eGType == ogr.wkbMultiPolygon or \
eGType == ogr.wkbGeometryCollection:
for i in range(poGeom.GetGeometryCount()):
SetZ(poGeom.GetGeometryRef(i), dfZ)
#**********************************************************************
# SetupTargetLayer()
#**********************************************************************
def SetZ (poGeom, dfZ ):
if poGeom is None:
return
eGType = wkbFlatten(poGeom.GetGeometryType())
if eGType == ogr.wkbPoint:
poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ)
elif eGType == ogr.wkbLineString or \
eGType == ogr.wkbLinearRing:
for i in range(poGeom.GetPointCount()):
poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ)
elif eGType == ogr.wkbPolygon or \
eGType == ogr.wkbMultiPoint or \
eGType == ogr.wkbMultiLineString or \
eGType == ogr.wkbMultiPolygon or \
eGType == ogr.wkbGeometryCollection:
for i in range(poGeom.GetGeometryCount()):
SetZ(poGeom.GetGeometryRef(i), dfZ)
#/************************************************************************/
#/* SetupTargetLayer() */
#/************************************************************************/
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 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 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 _get_vector_attributes(self,ds,attrs):
logging.debug("[FileMan][_get_vector_attributes]")
layer = ds.GetLayer()
# extent
extent = layer.GetExtent()
logging.debug("[FileMan][_get_vector_attributes] Extent: %s" % repr(extent) )
attrs["extent"] = (extent[0],extent[2],extent[1],extent[3])
# features count
attrs["features_count"] = layer.GetFeatureCount()
# geom type
ftype = layer.GetGeomType()
if ftype == ogr.wkbPolygon: #3
attrs["type"] = "polygon"
elif ftype == ogr.wkbPoint: #1
attrs["type"] = "point"
elif ftype == ogr.wkbLineString: #2
attrs["type"] = "line"
elif ftype == ogr.wkbPolygon25D:
attrs["type"] = "polygon"
else:
attrs["type"] = "none/unknown" # 0 or anything else
# srs
sr = layer.GetSpatialRef()
logging.debug("[FileMan][_get_vector_attributes] Spatial Reference from layer.GetSpatialRef() : %s" % repr(sr) )
if sr:
attrs["prj"]= self._get_prj(sr)
else:
attrs["prj"] = "unknown"
# Done
return attrs
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 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 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 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 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 createBoxFromLine(tmpGeom, ratio=1, halfWidth=-999, transformRequired=True, transform_WGS84_To_UTM='', transform_UTM_To_WGS84=''):
# create Polygon Box Oriented with the line
if transformRequired:
if transform_WGS84_To_UTM == '':
transform_WGS84_To_UTM, transform_UTM_To_WGS84 = createUTMTransform(tmpGeom)
tmpGeom.Transform(transform_WGS84_To_UTM)
# calculatuate Centroid
centroidX, centroidY, centroidZ = tmpGeom.Centroid().GetPoint()
lengthM = tmpGeom.Length()
if halfWidth ==-999:
halfWidth = lengthM/(2*ratio)
envelope=tmpGeom.GetPoints()
cX1 = envelope[0][0]
cY1 = envelope[0][1]
cX2 = envelope[1][0]
cY2 = envelope[1][1]
angRad = math.atan2(cY2-cY1,cX2-cX1)
d_X = math.cos(angRad-math.pi/2)*halfWidth
d_Y = math.sin(angRad-math.pi/2)*halfWidth
e_X = math.cos(angRad+math.pi/2)*halfWidth
e_Y = math.sin(angRad+math.pi/2)*halfWidth
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(cX1+d_X, cY1+d_Y)
ring.AddPoint(cX1+e_X, cY1+e_Y)
ring.AddPoint(cX2+e_X, cY2+e_Y)
ring.AddPoint(cX2+d_X, cY2+d_Y)
ring.AddPoint(cX1+d_X, cY1+d_Y)
polyGeom = ogr.Geometry(ogr.wkbPolygon)
polyGeom.AddGeometry(ring)
areaM = polyGeom.GetArea()
if transformRequired:
tmpGeom.Transform(transform_UTM_To_WGS84)
polyGeom.Transform(transform_UTM_To_WGS84)
return (polyGeom, areaM, angRad, lengthM)
def geoPolygonToPixelPolygonWKT(geom, inputRaster, targetSR, geomTransform, breakMultiPolygonGeo=True,
pixPrecision=2):
# Returns Pixel Coordinate List and GeoCoordinateList
polygonPixBufferList = []
polygonPixBufferWKTList = []
polygonGeoWKTList = []
if geom.GetGeometryName() == 'POLYGON':
polygonPix = ogr.Geometry(ogr.wkbPolygon)
for ring in geom:
# GetPoint returns a tuple not a Geometry
ringPix = ogr.Geometry(ogr.wkbLinearRing)
for pIdx in xrange(ring.GetPointCount()):
lon, lat, z = ring.GetPoint(pIdx)
xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)
xPix = round(xPix, pixPrecision)
yPix = round(yPix, pixPrecision)
ringPix.AddPoint(xPix, yPix)
polygonPix.AddGeometry(ringPix)
polygonPixBuffer = polygonPix.Buffer(0.0)
polygonPixBufferList.append([polygonPixBuffer, geom])
elif geom.GetGeometryName() == 'MULTIPOLYGON':
for poly in geom:
polygonPix = ogr.Geometry(ogr.wkbPolygon)
for ring in poly:
# GetPoint returns a tuple not a Geometry
ringPix = ogr.Geometry(ogr.wkbLinearRing)
for pIdx in xrange(ring.GetPointCount()):
lon, lat, z = ring.GetPoint(pIdx)
xPix, yPix = latlon2pixel(lat, lon, inputRaster, targetSR, geomTransform)
xPix = round(xPix, pixPrecision)
yPix = round(yPix, pixPrecision)
ringPix.AddPoint(xPix, yPix)
polygonPix.AddGeometry(ringPix)
polygonPixBuffer = polygonPix.Buffer(0.0)
if breakMultiPolygonGeo:
polygonPixBufferList.append([polygonPixBuffer, poly])
else:
polygonPixBufferList.append([polygonPixBuffer, geom])
for polygonTest in polygonPixBufferList:
if polygonTest[0].GetGeometryName() == 'POLYGON':
polygonPixBufferWKTList.append([polygonTest[0].ExportToWkt(), polygonTest[1].ExportToWkt()])
elif polygonTest[0].GetGeometryName() == 'MULTIPOLYGON':
for polygonTest2 in polygonTest[0]:
polygonPixBufferWKTList.append([polygonTest2.ExportToWkt(), polygonTest[1].ExportToWkt()])
return polygonPixBufferWKTList
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 DumpReadableGeometry( poGeometry, pszPrefix, options ):
if pszPrefix == None:
pszPrefix = ""
if 'DISPLAY_GEOMETRY' in options and EQUAL(options['DISPLAY_GEOMETRY'], 'SUMMARY'):
line = ("%s%s : " % (pszPrefix, poGeometry.GetGeometryName() ))
eType = poGeometry.GetGeometryType()
if eType == ogr.wkbLineString or eType == ogr.wkbLineString25D:
line = line + ("%d points" % poGeometry.GetPointCount())
print(line)
elif eType == ogr.wkbPolygon or eType == ogr.wkbPolygon25D:
nRings = poGeometry.GetGeometryCount()
if nRings == 0:
line = line + "empty"
else:
poRing = poGeometry.GetGeometryRef(0)
line = line + ("%d points" % poRing.GetPointCount())
if nRings > 1:
line = line + (", %d inner rings (" % (nRings - 1))
for ir in range(0,nRings-1):
if ir > 0:
line = line + ", "
poRing = poGeometry.GetGeometryRef(ir+1)
line = line + ("%d points" % poRing.GetPointCount())
line = line + ")"
print(line)
elif eType == ogr.wkbMultiPoint or \
eType == ogr.wkbMultiPoint25D or \
eType == ogr.wkbMultiLineString or \
eType == ogr.wkbMultiLineString25D or \
eType == ogr.wkbMultiPolygon or \
eType == ogr.wkbMultiPolygon25D or \
eType == ogr.wkbGeometryCollection or \
eType == ogr.wkbGeometryCollection25D:
line = line + "%d geometries:" % poGeometry.GetGeometryCount()
print(line)
for ig in range(poGeometry.GetGeometryCount()):
subgeom = poGeometry.GetGeometryRef(ig)
from sys import version_info
if version_info >= (3,0,0):
exec('print("", end=" ")')
else:
exec('print "", ')
DumpReadableGeometry( subgeom, pszPrefix, options)
else:
print(line)
elif 'DISPLAY_GEOMETRY' not in options or EQUAL(options['DISPLAY_GEOMETRY'], 'yes') \
or EQUAL(options['DISPLAY_GEOMETRY'], 'WKT'):
print("%s%s" % (pszPrefix, poGeometry.ExportToWkt() ))
return
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
poGeom = None
poDS = ogr.Open( pszDS, False )
if poDS is None:
return None
if pszSQL is not None:
poLyr = poDS.ExecuteSQL( pszSQL, None, None )
elif pszLyr is not None:
poLyr = poDS.GetLayerByName(pszLyr)
else:
poLyr = poDS.GetLayer(0)
if poLyr is None:
print("Failed to identify source layer from datasource.")
poDS.Destroy()
return None
if pszWhere is not None:
poLyr.SetAttributeFilter(pszWhere)
poFeat = poLyr.GetNextFeature()
while poFeat is not None:
poSrcGeom = poFeat.GetGeometryRef()
if poSrcGeom is not None:
eType = wkbFlatten(poSrcGeom.GetGeometryType())
if poGeom is None:
poGeom = ogr.Geometry( ogr.wkbMultiPolygon )
if eType == ogr.wkbPolygon:
poGeom.AddGeometry( poSrcGeom )
elif eType == ogr.wkbMultiPolygon:
for iGeom in range(poSrcGeom.GetGeometryCount()):
poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )
else:
print("ERROR: Geometry not of polygon type." )
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return None
poFeat = poLyr.GetNextFeature()
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return poGeom
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere):
poGeom = None
poDS = ogr.Open( pszDS, False )
if poDS is None:
return None
if pszSQL is not None:
poLyr = poDS.ExecuteSQL( pszSQL, None, None )
elif pszLyr is not None:
poLyr = poDS.GetLayerByName(pszLyr)
else:
poLyr = poDS.GetLayer(0)
if poLyr is None:
print("Failed to identify source layer from datasource.")
poDS.Destroy()
return None
if pszWhere is not None:
poLyr.SetAttributeFilter(pszWhere)
poFeat = poLyr.GetNextFeature()
while poFeat is not None:
poSrcGeom = poFeat.GetGeometryRef()
if poSrcGeom is not None:
eType = wkbFlatten(poSrcGeom.GetGeometryType())
if poGeom is None:
poGeom = ogr.Geometry( ogr.wkbMultiPolygon )
if eType == ogr.wkbPolygon:
poGeom.AddGeometry( poSrcGeom )
elif eType == ogr.wkbMultiPolygon:
for iGeom in range(poSrcGeom.GetGeometryCount()):
poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) )
else:
print("ERROR: Geometry not of polygon type." )
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return None
poFeat = poLyr.GetNextFeature()
if pszSQL is not None:
poDS.ReleaseResultSet( poLyr )
poDS.Destroy()
return poGeom
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 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
