Earth Map - error in running statistics due to complex polygon

  • Asked
  • Active
  • Viewed 2525 times
0

Dear Open Foris team,

I am trying to calculate the Land Productivity Dynamics with Earth Map using a complex polygon to be found in this folder: https://drive.google.com/drive/folders/1ITQTTHS2w4tVwZseIJurv0Ov5TNupQte The polygon has been generated following a river stream with a buffer of 12 km.

What do you suggest, shall I simplify it to perform the analytics?

0

Accepted answer

Problem:

The issue you're encountering while calculating analytics on Earth Map is not so much about the size of the area of your shapefile, but more about the complexity of the polygons themselves. Google Earth Engine, which Earth Map depends on, struggles with very complex polygons, particularly when there is clipping involved. This is compounded when the shapefile contains many small polygons and holes.

For example, the shapefile you're using includes many small islands in the rivers and holes (such as small lakes) outside the rivers.

Zoomed detail showing islands and holes:
Detail of islands and holes

These small polygons and holes cause performance issues during analysis. The best solution is to simplify the shapefile by removing them. You can do this using a tool like QGIS or with a simple Python script, which I'll guide you through below.


Solution: Simplify Your Shapefile with Python

Here’s how you can clean up your shapefile by removing small polygons (less than 100 hectares by default) and holes within polygons. This will make the shapefile easier to process on Earth Map and Google Earth Engine.

Steps:

  1. Install Python (if you don’t already have it).
  2. Run this command to install the necessary libraries:
   pip install geopandas
  1. Copy the code below into a file named preparePolygon.py.
  2. Run the script with:
   python preparePolygon.py

Important Notes: - Change the shapefile path to your file's location in line 29. - Check the CRS (Coordinate Reference System) for your area and adjust it in line 33. If you're unsure, feel free to ask on the forum. - You can modify the minimum polygon size (default is 100 hectares) in line 45. - The script outputs both a shapefile and a GeoJSON file. The GeoJSON can be used directly in Earth Map!

You can download the Python Script here


Python Script to Simplify Polygons:

import geopandas as gpd
from shapely.geometry import Polygon, MultiPolygon

# Function to remove holes from a polygon or multipolygon
def remove_holes(geometry):
    if isinstance(geometry, Polygon):
        return Polygon(geometry.exterior)  # Keep only the exterior (no holes)
    elif isinstance(geometry, MultiPolygon):
        return MultiPolygon([Polygon(p.exterior) for p in geometry.geoms])
    else:
        return geometry

# Function to round the coordinates of a polygon or multipolygon to a specific precision
def round_coordinates(geometry, precision=6):
    if isinstance(geometry, Polygon):
        exterior_rounded = [(round(x, precision), round(y, precision)) for x, y in geometry.exterior.coords]
        return Polygon(exterior_rounded)
    elif isinstance(geometry, MultiPolygon):
        return MultiPolygon([round_coordinates(p, precision) for p in geometry.geoms])
    else:
        return geometry

# Load your shapefile
shapefile_path = 'C:/Temp/Nigeria_StudyArea_EarthMap/Nigeria_StudyArea.shp'
gdf = gpd.read_file(shapefile_path)

# Reproject to UTM zone 31N (EPSG: 32631) for western Nigeria
gdf_projected = gdf.to_crs(epsg=32631)

# Remove holes from all geometries
gdf_projected['geometry'] = gdf_projected['geometry'].apply(remove_holes)

# Explode multipolygons into individual polygons for area filtering
gdf_exploded = gdf_projected.explode(index_parts=False)

# Calculate area in hectares (since UTM CRS is in meters)
gdf_exploded['area_ha'] = gdf_exploded.geometry.area / 10000  # Convert square meters to hectares

# Filter out polygons smaller than 100 hectares
gdf_filtered = gdf_exploded[gdf_exploded['area_ha'] >= 100]

# Group polygons back into multipolygons where necessary
gdf_final = gdf_filtered.dissolve(by=gdf_filtered.index)

# Reproject the final GeoDataFrame to EPSG:4326 (WGS84)
gdf_final = gdf_final.to_crs(epsg=4326)

# Round the geometry coordinates to 6 decimal places
gdf_final['geometry'] = gdf_final['geometry'].apply(lambda geom: round_coordinates(geom, precision=6))

# Save the filtered shapefile
shapefile_output = 'shapefileWithNoHolesOrSmallPolygons.shp'
gdf_final.to_file(shapefile_output)

# Save the filtered GeoDataFrame as a GeoJSON
geojson_output = 'shapefileWithNoHolesOrSmallPolygons.geojson'
gdf_final.to_file(geojson_output, driver='GeoJSON')

print(f"Filtered shapefile saved to: {shapefile_output}")
print(f"Filtered GeoJSON saved to: {geojson_output}")

Output Files:

  • Shapefile: shapefileWithNoHolesOrSmallPolygons.shp
  • GeoJSON: shapefileWithNoHolesOrSmallPolygons.geojson

You can now upload the cleaned GeoJSON directly into Earth Map for better performance.

Let us know if this solution works for you! We’ve tested it on your shapefile, and here are the results:

Result

Your answer

Markdown text. Shortcuts: Ctrl+B bold, Ctrl+I italic, Ctrl+K link, Ctrl+E code, Ctrl+Shift+V paste as plain text. Tab indents list items and code; press Escape, then Tab, to leave the text box.

Insert an image

Say what the image shows, for people who can't see it.

Keyboard shortcuts

Anywhere

?
Show this list
/
Search
g h
Go to the home page
g q
Go to the questions
g t
Go to the tags
g u
Go to the users
a
Ask a question

Lists of questions

j k
Next or previous question
Enter
Open the question

Question pages

j k
Next or previous post

Shortcuts never apply while you type in a text field. Members can turn them off in their preferences.