TransWikia.com

Is it possible to create a shapefile from a geometry with Python?

Geographic Information Systems Asked on August 24, 2021

I use a GEE asset to describe the geometry of a country. I can transform this geometry into different format such as shapely.geometry or GeoJSON:

asset = 'users/bornToBeAlive/aoi_PU'

aoi = ee.FeatureCollection(asset)
aoiJson = geemap.ee_to_geojson(aoi)
aoiShp = sg.shape(aoiJson['features'][0]['geometry'])

rootDir = os.path.expanduser('~') + '/'

Unfortunately the gdal.warp function requires a shapefile as input. Is there a way to export these structure to a shapefile ?

EDIT

# Now convert it to a shapefile with OGR    
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(rootDir + 'my.shp')
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()

# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', 123)

# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(aoiShp.wkb)
feat.SetGeometry(geom)

layer.CreateFeature(feat)
feat = geom = None  # destroy these

# Save and close everything
ds = layer = feat = geom = None

I tried to use this solution but I obtain a very small version of my geometry centered in [0,0] when I open it in QGIS, am I missing something ?

One Answer

So, in order to create a shapefile from an ee asset I propose the following solution, feel free to improve it :

rootDir = os.path.expanduser('~') + '/'

asset = 'users/bornToBeAlive/aoi_PU'

aoi = ee.FeatureCollection(asset)
aoiJson = geemap.ee_to_geojson(aoi)
aoiShp = sg.shape(aoiJson['features'][0]['geometry'])

Then you need to create your shapefile (it will be encoded in UTF-8)

# Now convert it to a shapefile with OGR    
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource(rootDir + 'my.shp')
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()

# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', 123)

# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(aoiShp.wkb)
feat.SetGeometry(geom)

layer.CreateFeature(feat)
feat = geom = None  # destroy these

# Save and close everything
ds = layer = feat = geom = None

Finally add a projection file to avoid projection error when displaying or using:

#add the spatial referecence
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)

spatialRef.MorphToESRI()
file = open(rootDir + 'my.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

and then you are good to go !

Correct answer by Pierrick Rambaud on August 24, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP