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类GetDriverByName()的实例源码
def parsed(args):
kml_file=args.geo
def kml2geojson(kml_file):
drv = ogr.GetDriverByName('KML')
kml_ds = drv.Open(kml_file)
for kml_lyr in kml_ds:
for feat in kml_lyr:
outfile=feat.ExportToJson()
geom2=str(outfile).replace(", 0.0",'')
with open(args.loc+'./kmlout.geojson','w') as csvfile:
writer=csv.writer(csvfile)
writer.writerow([geom2])
kml2geojson(args.geo)
raw= open(args.loc+'./kmlout.geojson')
for line in raw:
fields=line.strip().split(":")[3]
f2=fields.strip().split("}")[0]
filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
with open(args.loc+'./aoi.json', 'w') as outfile:
outfile.write(filenames)
outfile.close()
def parsed(args):
kml_file=args.geo
def kml2geojson(kml_file):
drv = ogr.GetDriverByName('KML')
kml_ds = drv.Open(kml_file)
for kml_lyr in kml_ds:
for feat in kml_lyr:
outfile=feat.ExportToJson()
geom2=str(outfile).replace(", 0.0",'')
with open(args.loc+'./kmlout.geojson','w') as csvfile:
writer=csv.writer(csvfile)
writer.writerow([geom2])
kml2geojson(args.geo)
raw= open(args.loc+'./kmlout.geojson')
for line in raw:
fields=line.strip().split(":")[3]
f2=fields.strip().split("}")[0]
filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
with open(args.loc+'./aoi.json', 'w') as outfile:
outfile.write(filenames)
outfile.close()
def import_summary_geojson(geojsonfilename, removeNoBuildings=True):
# driver = ogr.GetDriverByName('geojson')
datasource = ogr.Open(geojsonfilename, 0)
layer = datasource.GetLayer()
print(layer.GetFeatureCount())
buildingList = []
for idx, feature in enumerate(layer):
poly = feature.GetGeometryRef()
if poly:
if removeNoBuildings:
if feature.GetField('BuildingId') != -1:
buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
'poly': feature.GetGeometryRef().Clone()})
else:
buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
'poly': feature.GetGeometryRef().Clone()})
return buildingList
def import_chip_geojson(geojsonfilename, ImageId=''):
# driver = ogr.GetDriverByName('geojson')
datasource = ogr.Open(geojsonfilename, 0)
if ImageId=='':
ImageId = geojsonfilename
layer = datasource.GetLayer()
print(layer.GetFeatureCount())
polys = []
BuildingId = 0
for idx, feature in enumerate(layer):
poly = feature.GetGeometryRef()
if poly:
BuildingId = BuildingId + 1
polys.append({'ImageId': ImageId, 'BuildingId': BuildingId,
'poly': feature.GetGeometryRef().Clone()})
return polys
def createRasterFromGeoJson(srcGeoJson, srcRasterFileName, outRasterFileName):
NoData_value = 0
source_ds = ogr.Open(srcGeoJson)
source_layer = source_ds.GetLayer()
srcRaster = gdal.Open(srcRasterFileName)
# Create the destination data source
target_ds = gdal.GetDriverByName('GTiff').Create(outRasterFileName, srcRaster.RasterXSize, srcRaster.RasterYSize, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform(srcRaster.GetGeoTransform())
target_ds.SetProjection(srcRaster.GetProjection())
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)
# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
band.FlushCache()
def parsed(args):
kml_file=args.geo
def kml2geojson(kml_file):
drv = ogr.GetDriverByName('KML')
kml_ds = drv.Open(kml_file)
for kml_lyr in kml_ds:
for feat in kml_lyr:
outfile=feat.ExportToJson()
geom2=str(outfile).replace(", 0.0",'')
with open(args.loc+'./kmlout.geojson','w') as csvfile:
writer=csv.writer(csvfile)
writer.writerow([geom2])
kml2geojson(args.geo)
raw= open(args.loc+'./kmlout.geojson')
for line in raw:
fields=line.strip().split(":")[3]
f2=fields.strip().split("}")[0]
filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
with open(args.loc+'./aoi.json', 'w') as outfile:
outfile.write(filenames)
outfile.close()
def mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32):
"""Create a new GDAL Dataset in memory
Useful for various applications that require a Dataset
"""
#These round down to int
#dst_ns = int((extent[2] - extent[0])/res)
#dst_nl = int((extent[3] - extent[1])/res)
#This should pad by 1 pixel, but not if extent and res were calculated together to give whole int
dst_ns = int((extent[2] - extent[0])/res + 0.99)
dst_nl = int((extent[3] - extent[1])/res + 0.99)
m_ds = gdal.GetDriverByName('MEM').Create('', dst_ns, dst_nl, 1, dtype)
m_gt = [extent[0], res, 0, extent[3], 0, -res]
m_ds.SetGeoTransform(m_gt)
if srs is not None:
m_ds.SetProjection(srs.ExportToWkt())
return m_ds
#Modify proj/gt of dst_fn in place
def __init__(self, filename=None, driver="ESRI Shapefile"):
if driver not in ["ESRI Shapefile", "Memory"]:
raise IOError("driver not supported")
if filename is None:
driver = "Memory"
else:
self.filename = filename
self.driver = ogr.GetDriverByName(driver)
self.vector = self.driver.CreateDataSource("out") if driver == "Memory" else self.driver.Open(filename)
nlayers = self.vector.GetLayerCount()
if nlayers > 1:
raise IOError("multiple layers are currently not supported")
elif nlayers == 1:
self.init_layer()
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 generate_tiles(self):
"""INCOMPLETE Split the calculated WPbm in 100 tiles facilitating the export
Returns:
TYPE: Description
"""
driver = ogr.GetDriverByName('ESRI Shapefile')
dir_shps = "tiles"
os.chdir(dir_shps)
file_shps = glob.glob("*.shp")
allExportWPbm = []
file_names = []
for file_shp in file_shps:
dataSource = driver.Open(file_shp, 0)
if dataSource is None:
sys.exit(('Could not open {0}.'.format(file_shp)))
else:
layer = dataSource.GetLayer(0)
extent = layer.GetExtent()
active_file = "tile_" + str(file_shp.split('.')[0]).split("_")[3]
file_names.append(active_file)
low_sx = extent[0], extent[3]
up_sx = extent[0], extent[2]
up_dx = extent[1], extent[2]
low_dx = extent[1], extent[3]
cut = [list(low_sx), list(up_sx), list(up_dx), list(low_dx)]
Export_WPbm = {
"crs": "EPSG:4326",
"scale": 250,
'region': cut}
allExportWPbm.append(Export_WPbm)
return allExportWPbm, file_names
def kml2geojson(kml_file):
'''
From Gist.
Transform kml (from Google Maps) to the geojson format.
'''
drv = ogr.GetDriverByName('KML')
kml_ds = drv.Open(kml_file)
for kml_lyr in kml_ds:
for feat in kml_lyr:
print feat.ExportToJson()
def get_bbox(shapefile):
""" Gets the bounding box of a shapefile in EPSG 4326.
If shapefile is not in WGS84, bounds are reprojected.
"""
driver = ogr.GetDriverByName('ESRI Shapefile')
# 0 means read, 1 means write
data_source = driver.Open(shapefile, 0)
# Check to see if shapefile is found.
if data_source is None:
failure('Could not open %s' % (shapefile))
else:
layer = data_source.GetLayer()
shape_bbox = layer.GetExtent()
spatialRef = layer.GetSpatialRef()
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
# this check for non-WGS84 projections gets some false positives, but that's okay
if target.ExportToProj4() != spatialRef.ExportToProj4():
transform = osr.CoordinateTransformation(spatialRef, target)
point1 = ogr.Geometry(ogr.wkbPoint)
point1.AddPoint(shape_bbox[0], shape_bbox[2])
point2 = ogr.Geometry(ogr.wkbPoint)
point2.AddPoint(shape_bbox[1], shape_bbox[3])
point1.Transform(transform)
point2.Transform(transform)
shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]
return shape_bbox
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 _get_ft_ds():
refresh_token = OAuth2().get_refresh_token()
ft_driver = ogr.GetDriverByName('GFT')
return ft_driver.Open('GFT:refresh=' + refresh_token, True)
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 __init__(self, country_file):
driver = ogr.GetDriverByName('ESRI Shapefile')
self.countryFile = driver.Open(country_file)
self.layer = self.countryFile.GetLayer()
def open_vector(filename, driver=None):
"""Open vector file, return gdal.Dataset and OGR.Layer
.. warning:: dataset and layer have to live in the same context,
if dataset is deleted all layer references will get lost
.. versionadded:: 0.12.0
Parameters
----------
filename : string
vector file name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
layer : ogr.Layer
layer
"""
dataset = gdal.OpenEx(filename)
if driver:
gdal.GetDriverByName(driver)
layer = dataset.GetLayer()
return dataset, layer
def open_shape(filename, driver=None):
"""Open shapefile, return gdal.Dataset and OGR.Layer
.. warning:: dataset and layer have to live in the same context,
if dataset is deleted all layer references will get lost
.. versionadded:: 0.6.0
Parameters
----------
filename : string
shapefile name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
layer : ogr.Layer
layer
"""
if driver is None:
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(filename)
if dataset is None:
print('Could not open file')
raise IOError
layer = dataset.GetLayer()
return dataset, layer
def open_raster(filename, driver=None):
"""Open raster file, return gdal.Dataset
.. versionadded:: 0.6.0
Parameters
----------
filename : string
raster file name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
"""
dataset = gdal.OpenEx(filename)
if driver:
gdal.GetDriverByName(driver)
return dataset
def read_safnwc(filename):
"""Read MSG SAFNWC hdf5 file into a gdal georeferenced object
Parameters
----------
filename : string
satellite file name
Returns
-------
ds : gdal.DataSet
with satellite data
"""
root = gdal.Open(filename)
ds1 = gdal.Open('HDF5:' + filename + '://CT')
ds = gdal.GetDriverByName('MEM').CreateCopy('out', ds1, 0)
# name = os.path.basename(filename)[7:11]
try:
proj = osr.SpatialReference()
proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
except Exception:
raise NameError("No metadata for satellite file %s" % filename)
geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
ds.SetProjection(proj.ExportToWkt())
ds.SetGeoTransform([float(x) for x in geotransform])
return ds
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 _getOgrTestLayer(self):
dest = self._copyTestLayer()
driver = ogr.GetDriverByName('GPKG')
dataSource = driver.Open(dest, 1)
return dataSource
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)
