InsideDarkWeb.com

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

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 ?

Geographic Information Systems Asked on November 21, 2021

1 Answers

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 !

Answered by Pierrick Rambaud on November 21, 2021

Add your own answers!

Related Questions

How to open shapefile dataset using ArcGIS Pro SDK .NET?

1  Asked on June 24, 2021 by pavel-butusov

   

GRASS 6.4.3 Crashes on Startup; libintl-8.dll cited

4  Asked on June 24, 2021 by darin

     

Plotting polygons as separate plots using Python

3  Asked on June 23, 2021 by ralt_jol

       

Seeking formula for Distance Calculation on Mars

0  Asked on June 23, 2021 by douglas-marsh

     

Ask a Question

Get help from others!

© 2021 InsideDarkWeb.com. All rights reserved.