I have GeoTIFF images from https://esa-worldcover.org/en

My goal is to get a list of latitude / longitude / value from that image. The value is a numeric attribute which represents a class as explained in https://worldcover2020.esa.int/data/docs/WorldCover_PUM_V1.1.pdf

I could get the value from each pixel by doing:

import rasterio
ds = rasterio.open(img)

However, when I'm trying to link those values from lat and lon (after transforming x and y values to lat and lon), the shape doesn't fit.

The code I've used for transforming x and y into coordinates is:

import rasterio
import rasterio.features
import rasterio.warp

with rasterio.open('ESA_WorldCover_10m_2020_v100_N39E000_Map.tif') as dataset:

    # Read the dataset's valid data mask as a ndarray.
    mask = dataset.dataset_mask()
    # Extract feature shapes and values from the array.
    for geom, val in rasterio.features.shapes(
            mask, transform=dataset.transform):
        # Transform shapes from the dataset's own coordinate
        # reference system to CRS84 (EPSG:4326).
        geom = rasterio.warp.transform_geom(
            dataset.crs, 'EPSG:4326', geom, precision=8) 


I only get 13213 values. How is it possible? The map can be downloaded in here

  • rasterio.features.shapes returns connected groups of pixels of the same value, not individual pixels. Just convert your tiff to XYZ e.g gdal_translate -of XYZ input.tif output.xyz then read it in with pandas.
    – user2856
    Apr 12, 2022 at 23:44

1 Answer 1


This should do the job...

(There's no need to obtain all coordinates of the raster first (which would be more than 1 billion in this case...)

import rasterio
dat = rasterio.open(r"ESA_WorldCover_10m_2020_v100_N45W120_Map.tif")
# read all the data from the first band
z = dat.read()[0]

# check the crs of the data
# >>> CRS.from_epsg(4326)

# check the bounding-box of the data
# >>> Out[49]: BoundingBox(left=-120.0, bottom=45.0, right=-117.0, top=48.0)

# since the raster is in regular lon/lat grid (4326) we can use 
# `dat.index()` to identify the index of a given lon/lat pair
# (e.g. it expects coordinates in the native crs of the data)

def getval(lon, lat):
    idx = dat.index(lon, lat, precision=1E-6)    
    return dat.xy(*idx), z[idx]

getval(-118, 46)
# >>> ((-117.99995833333334, 46.00004166666667), 10)

... and so we found that (-118, 46) seems to be covered by Trees :-)

enter image description here

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.