1

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 ?

6
  • does this work? gis.stackexchange.com/questions/52705/…
    – risail
    Commented Jul 17, 2020 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 ? Commented Jul 22, 2020 at 7:59
  • What projection (EPSG) do you use?
    – Zoltan
    Commented Jul 22, 2020 at 8:09
  • To display it in QGIS I need to manually add 'ESPG: 4326' is it possible to add it to ds ? Commented Jul 22, 2020 at 8:12
  • gdal.warp function requires a shapefile as input is the function the same as gdalwarp?
    – nmtoken
    Commented Jul 22, 2020 at 8:24

1 Answer 1

2

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 !

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.