Geographic Information Systems Asked on November 21, 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 ?
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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP