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

Question

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 ?

does this work? https://gis.stackexchange.com/questions/52705/how-to-write-shapely-geometries-to-shapefiles — risail, Jul 17 '20 at 13:10
I've tested the method of this link but when I open my sapefile in QGIS to check it's the appropriate geometry but super small and centered at the [0,0] point, did I miss something ? — Pierrick Rambaud, Jul 22 '20 at 07:59
To display it in QGIS I need to manually add 'ESPG: 4326' is it possible to add it to ds ? — Pierrick Rambaud, Jul 22 '20 at 08:12
gdal.warp function requires a shapefile as input is the function the same as gdalwarp? — nmtoken, Jul 22 '20 at 08:24

score 2 · Accepted Answer · answered Jul 22 '20 at 08:26

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 !

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

Create a new feature (attribute and geometry)

Make a geometry, from Shapely object

Save and close everything

1 Answers1

Create a new feature (attribute and geometry)

Make a geometry, from Shapely object

Save and close everything