ESA WorldCereal - Getting Locations for Global Field-Level Monitoring in Python


With the release of the ESA WorldCereal system, we now have (for the first time ever) the precise location of essentially every winter cereal, spring cereal and maize farm in the world.

To put the precision in perspective, the best datasets available to date were at 10km and ~5km resolutions via SPAM2010 and GEOGLAM, respectively. WorldCereal is 10m resolution – that’s 1000 times more precise. In fact, this is even better than the 30m resolution that the USDA has on offer with the USDA CDL (Cropland Data Layer).

WorldCereal is monumental because it is the missing link that can enable global monitoring of the entire world cereal situation. That said, the key to unlocking that potential is to transform WorldCereals from its cumbersome native format (.tif / geotiff rasters) to point samples in Parquet or GeoParquet.

Once that is done, there is no limit to what data scientists can do with this treasure trove. For example, measuring crop health for any (or all!) of the world’s cereal farms is imminently doable. (Use NDVI as your metric and you can do this with our API in a single Python program.)

To that end then, this tutorial shows you exactly how to convert WorldCereals into points samples in Parquet or GeoParquet. It will show how to:

  • Generate any number of global maize, winter cereals, and spring cereals farm point locations (longitude, latitude coordinate pairs) from the .tif classification mask files
  • Export the locations to cloud-native (Parquet & GeoParquet) and standard tabular (CSV) formats

Parquet and GeoParquet enable maximum interoperability with data science workflows and the modern stack, including cloud storage like AWS, GCP, and Azure, cloud data warehouses like Snowflake, BigQuery, RedShift, and DataBricks, and data science / machine learning platforms like pandas, numpy, scikit-learn, and more.

About WorldCereal

The ESA WorldCereal dataset consists of 10m resolution global crop type maps (classification masks) covering maize, winter cereals and spring cereals. For maize and winter cereals, there are 106 .tif (raster) files corresponding to the 106 Agricultural Ecological Zones (AEZs) across the globe. Spring cereals are only covered for 21 of the 106 AEZs. Find the article (Van Tricht et al., 2023, in pre-print) describing the full methodology here. Note, while not covered in this tutorial, the WorldCereal dataset also has classification mask products for seasonal irrigated cropland and temporary crop extent which can be processed in the same manner as the 3 covered here.

Note: this tutorial is also available as a single python notebook. You can download it on Github below:

Preconstructed Data

If you don’t want to generate your own locations by running this tutorial, you can download the output at the following links:

The output datasets consist of 6,827,917 generated point locations and have the following information for each: Crop, Country (Admin 0), State (Admin 1), Region (World Bank), Subregion, Continent, and aez_id. See below: 

Requisite Data

To begin you need to pull the WorldCereals datasets:

Download all of the above, unzip and move the contents to a folder called ‘data’ in your working directory. (The full WorldCereals dataset is available here. For smaller, user-specified regions it can also be accessed through the OpenEO API.) 

Requisite Modules

There are a number of libraries that need to be imported. The first group is required for generating the point locations from the raster geotiff files (.tif). The second group is required for tagging the point locations with various administrative parameters. 

# Required for location generation
import os
import geopandas as gpd
import rasterio
import numpy as np
import pandas as pd

# Required for location tagging
import cartopy
import as ccrs
import cartopy.feature as cf
import as shpreader
from shapely.geometry import Point, Polygon

Define Functions

The .tif raster files are very large so we first define a function to read them in chunks to avoid any memory-related issues.

def read_raster_in_chunks(dataset):
for ji, window in dataset.block_windows(1):
band_window =, indexes=1)
yield band_window, ji

Demo vs Full Run

Next we introduce a boolean constant called “demo_mode”. In demo-mode, running the program will generate Maize farm locations for just 3 of the 106 AEZs (covering Mexico, Guatemala, El Salvador, Honduras, and North-Western Nicaragua). It should take  less than 10min to run for these AEZs. If you wish to run the code on different AEZs covering different regions than specified, modify the aez_id_list. 

With demo-mode set to False, running the program will generate farm locations for maize, winter cereals, and spring cereals for all 106 AEZs (where available). Note running the notebook in full-mode can take tens of hours (probably between 10 and 20).

# Specify if the run is a demo
demo_mode = True # Change to False for a full run

Now we can proceed with point generation. The code below will iterate through each of the raster .tif files with rasterio for the list of specified AEZs, reading and processing each the .tif files in chunks. For each chunk, we calculate the total number of pixels classified as being a specific crop and then randomly select 0.01% of these pixels in which we will randomly generate points. This method allows us to generate points that will be scaled by the density of the crops in the original .tif files. Note that the outputs we show in this tutorial were produced with demo-mode set to False.

# Specify if the run is a demo
demo_mode = True # Change to False for a full run

# Read the geojson containing the aez footprints as polygons
aez = gpd.read_file('data/WorldCereal_AEZ.geojson')

# List of all aez_ids
aez_ids = aez['aez_id'].unique()

# Create a list to hold dataframes
dfs = []

# Define the directories and crop_types based on whether it is a demo run or not
if demo_mode:
    directories = ['data/WorldCereal_2021_tc-maize-main_maize_classification']
    crop_types = {'data/WorldCereal_2021_tc-maize-main_maize_classification': 'Maize'}
    directories = [
    crop_types = {
        'data/WorldCereal_2021_tc-maize-main_maize_classification': 'Maize',
        'data/WorldCereal_2021_tc-wintercereals_wintercereals_classification': 'Winter Cereals',
        'data/WorldCereal_2021_tc-springcereals_springcereals_classification': 'Spring Cereals'

# Determine the aez_id list for processing
aez_id_list = [15081,19093,33117] if demo_mode else aez_ids

# Iterate over all directories
for directory in directories:

    print(f"Processing: {crop_types[directory]} ({directory})")

    # Iterate over all unique aez_ids
    for aez_id in aez_id_list:

        print(f"Processing aez_id: {aez_id}")

        # Get list of all .tif files in the directory
        all_files = os.listdir(directory)

        for file in all_files:
                    # Check if the file starts with the aez_id and ends with .tif
                    if file.startswith(str(aez_id)) and file.endswith('.tif'):
                        # Full path to the file
                        filename = os.path.join(directory, file)

                        # Open the raster
                        dataset =

                        # Read the band in chunks
                        chunks = read_raster_in_chunks(dataset)

                        # Process one chunk at a time
                        for i, (band_chunk, ji) in enumerate(chunks):

                            # Find the indices where band_chunk has values of 100
                            indices = np.where(band_chunk == 100)

                            # Determine the total number of classified pixels
                            total_classified = len(indices[0])

                            # Determine the proportion of classified pixels to sample
                            proportion = 0.0001  # For example, to sample 0.01% of classified pixels

                            # Determine the number of samples to draw
                            num_samples = int(np.round(total_classified * proportion))

                            # Skip this chunk if there are fewer than 1 pixels with a value of 100
                            if num_samples < 1:

                            # Randomly sample num_samples points from the classified indices 
                            sample_indices = np.random.choice(total_classified, size=num_samples, replace=False)

                            # Get the row and column indices
                            rows = indices[0][sample_indices] + ji[0]*band_chunk.shape[0]
                            cols = indices[1][sample_indices] + ji[1]*band_chunk.shape[1]

                            # Transform the row and column indices to the original CRS (EPSG:4326)
                            x_coords, y_coords = dataset.xy(rows, cols)

                            data = {'lon': x_coords, 'lat': y_coords, 'aez_id': aez_id, 'crop': crop_types[directory]}

                            df = pd.DataFrame(data)

                            # Append the df to dfs

# Concatenate all the dataframes in dfs
final_df = pd.concat(dfs, ignore_index=True)

Executing the above generates 7M points (for all 3 crop types across all applicable AEZs):

Tag Administrative Regions

The next step is to tag each of the points by the administrative boundaries they fall inside. This makes the dataset versatile for multiple use cases.

We use the Natural Earth Features country and state-level administrative datasets to tag the farm locations generated above. For simplicity, we will access the datasets through Cartopy, but the full list of datasets may also be accessed here. Note that you can of course tag locations using any boundary file that is appropriate for your use case. 

# shpfilename = the path to your shapefile
ne_level_0 = shpreader.natural_earth(resolution='10m', category='cultural', name='admin_0_countries')
ne_level_1 = shpreader.natural_earth(resolution='10m', category='cultural', name='admin_1_states_provinces')

# Read the shapefile using GeoPandas
countries_gdf = gpd.read_file(ne_level_0)
states_gdf = gpd.read_file(ne_level_1)

Countries (Admin Level 0), Regions, Subregions, Continents

# Convert the df DataFrame to a GeoDataFrame
geometry = [Point(xy) for xy in zip(final_df['lon'], final_df['lat'])]
merged_geo = gpd.GeoDataFrame(final_df, geometry=geometry)

# Set the CRS of the merged_geo GeoDataFrame to EPSG:4326 (WGS 84) = "EPSG:4326"

# Transform the CRS of merged_geo to match the CRS of county_boundaries
merged_geo = merged_geo.to_crs(

# Perform a spatial join
joined = gpd.sjoin(merged_geo, countries_gdf, predicate='within', how='left')

States (Admin Level 1)

# Convert the df DataFrame to a GeoDataFrame
geometry = [Point(xy) for xy in zip(final_df['lon'], final_df['lat'])]
merged_geo = gpd.GeoDataFrame(final_df, geometry=geometry)

# Set the CRS of the merged_geo GeoDataFrame to EPSG:4326 (WGS 84) = "EPSG:4326"

# Transform the CRS of merged_geo to match the CRS of county_boundaries
merged_geo = merged_geo.to_crs(

# Perform a spatial join
joined = gpd.sjoin(merged_geo, states_gdf, predicate='within', how='left')

# Remove duplicate index after spatial join
joined = joined.loc[~joined.index.duplicated(keep='first')]

# # Copy the required columns back to the df DataFrame
final_df['State (Admin 1)'] = joined['name']

The DataFrame which originally contained only points locations and the corresponding crops is now tagged for the various administrative boundaries. Inspecting it visually:

Remove excess points generated in regions with AEZ overlap 

Because AEZ regions sometimes overlap each other, our random sampling methodology will yield excess points in areas of overlap. This can actually be seen visually in the image above. The density of points in these areas of overlap is around 40% higher than what it should be. There is an easy fix for this: delete 40% of the points that fall in regions where AEZ overlap occurs. To do this, use the WorldCereal_AEZ.geojson file containing the AEZ footprint polygons. 

# Convert the DataFrame to a GeoDataFrame
geometry = [Point(xy) for xy in zip(final_df['lon'], final_df['lat'])]
geo_df = gpd.GeoDataFrame(final_df, geometry=geometry)

# Set CRS for geo_df to match aez
geo_df.set_crs(epsg=4326, inplace=True)

# Read the geojson to a GeoDataFrame
aez = gpd.read_file('data/WorldCereal_AEZ.geojson')

# Create a spatial join between the points and the polygons to determine which points fall in overlapping polygons
overlaps = gpd.sjoin(geo_df, aez, how='inner', predicate='within')

# Group points by 'aez_id' and find the groups that have more than one polygon (overlapping areas)
overlapping_points = overlaps.groupby(['lon', 'lat']).filter(lambda x: len(x) > 1)

# Sample 40% of the rows 
overlapping_points_sampled = overlapping_points.sample(frac=0.40)

# If you want to reset index
overlapping_points_sampled.reset_index(drop=True, inplace=True)

# Create a key in both dataframes for merging
overlapping_points_sampled['key'] = overlapping_points_sampled['lon'].astype(str) + '-' + overlapping_points_sampled['lat'].astype(str)
final_df['key'] = final_df['lon'].astype(str) + '-' + final_df['lat'].astype(str)

# Perform a left merge to identify rows in final_df that are not present in overlapping_points_sampled
merged_df = pd.merge(final_df, overlapping_points_sampled[['key']], on='key', how='left', indicator=True)

# Filter rows that are only present in final_df (not present in overlapping_points_sampled)
final_df_filtered = merged_df[merged_df['_merge'] == 'left_only']

# Remove the key and _merge columns
final_df_filtered = final_df_filtered.drop(columns=['key', '_merge'])

final_df_filtered = final_df_filtered[final_df_filtered['aez_id'] != 2008]

By removing the excess points we have reduced the size of output by nearly 400,000 points:


Upon visual inspection, we can see that the excess points that were generated in areas of AEZ overlap have been eliminated: 

Filter Out Problem Areas

During the generation process, points can show up along the edges of some of the AEZs in regions that do not correspond to farms. We manually drew rectangular polygons around all areas where these points occurred, and saved the polygons in geojson format in a list.

problem_areas = [{"type":"Polygon","coordinates":[[[-91.41645262581089,54.84799699573268],[-91.41645262581089,54.40865331122363],[-88.50502911906254,54.40865331122363],[-88.50502911906254,54.84799699573268],[-91.41645262581089,54.84799699573268]]]},

We use the list of polygons to filter out any points that may have been generated within these known problem areas.

# Create the list of Polygon objects from the geojson data
problem_polygons = [Polygon(area['coordinates'][0]) for area in problem_areas]

# Function to check if a point is in any of the problem areas
def is_in_problem_area(row):
    point = Point(row['lon'], row['lat'])
    return any(polygon.contains(point) for polygon in problem_polygons)

# Apply the function to each row in the dataframe
final_df_filtered['is_in_problem_area'] = final_df_filtered.apply(is_in_problem_area, axis=1)

# Filter the dataframe to remove points in problem areas
final_df_filtered = final_df_filtered[final_df_filtered['is_in_problem_area'] == False]

# You can drop the 'is_in_problem_area' column afterwards if you wish
final_df_filtered = final_df_filtered.drop(columns=['is_in_problem_area'])

final_df_filtered = final_df_filtered.reset_index()
final_df_filtered = final_df_filtered.drop(columns=['index'])


Final Processing

We’ll do some final formatting to the dataset we’ve generated. Specifically, we reorder the columns into something that makes a little bit more sense and capitalize the ‘crop’ column.

final_df_filtered = final_df_filtered[['lon', 'lat', 'crop', 'Country (Admin 0)', 'State (Admin 1)', 'Region (World Bank)', 'Subregion', 'Continent', 'aez_id']]
final_df_filtered = final_df_filtered.rename(columns={'crop': 'Crop'})


Saving Final Locations Dataset

Finally, we will save the locations dataset in both cloud-native (Parquet & GeoParquet) and standard tabular (CSV) formats.

Parquet (.parquet)

final_df_filtered.to_parquet('world_cereals_scaled_final_all.parquet') # 101 MB

GeoParquet (.geoparquet)

# Creating the geometry column from the coordinates
geometry = [Point(xy) for xy in zip(final_df_filtered['lon'], final_df_filtered['lat'])]

# Creating a GeoDataFrame by adding the geometry column to the original DataFrame
final_gdf = gpd.GeoDataFrame(final_df_filtered, geometry=geometry)

final_gdf.drop(columns=['lon','lat']).to_parquet('world_cereals_scaled_final_all.geoparquet') # 100 MB

CSV (.csv)

For CSV, we will export both the full dataset we generated, and a random 3M subset of the points.

# final_df_filtered.sample(n=3000000, random_state=1).to_csv('world_cereals_scaled_final_3M.csv')

The ESA WorldCereal dataset will have massive implications for our ability to monitor crops no matter where they are in the world. With WorldCereal as point locations in the cloud-native Parquet and GeoParquet formats, these datasets are immediately unlocked for data science workflows and those wishing to do time series analytics with remote sensing data. From building more accurate crop yield models than has ever been possible, to constructing higher-resolution crop calendars globally, point sampling WorldCereals opens up all kinds of analysis possibilities. .