Using Rasterio to index by traditional longitude and latitude with Molleweide Projection GeoTIFF

Question

I'm currently trying to read carbon footprint information from a GeoTIFF, I've never worked with geographic data before.

The link is http://worldmrio.com/citiesmaps/GGMCF_v1.0.tif

I was wondering how I would index the nd-array with longitude and latitude values that I've scraped for some cities. What I want is the actual value at each longitude and latitude.

Additionally, would there be a problem with the projection format?

print(dataset.crs)

gives me

PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS    84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
{'type': 'Polygon', 'coordinates': [[(-18040095.3812882, 9020047.84316439), (-18040095.3812882, -9020047.84316439), (18040095.3812882, -9020047.84316439), (18040095.3812882, 9020047.84316439), (-18040095.3812882, 9020047.84316439)]]}

user2856 · Accepted Answer

You'll need to reproject your lon, lats to the Mollweide CRS (or reproject your raster to EPSG:4326 which may not be the best option).

There's numerous ways to reproject your coords, (geopandas, pyproj) the example below uses fiona.transform.

To extract the raster values at your coordinates, you can use rasterio.sample.

For example:

import fiona.transform
import rasterio.sample

def reproject_coords(src_crs, dst_crs, coords):
    xs = [c[0] for c in coords]
    ys = [c[1] for c in coords]
    xs, ys = fiona.transform.transform(src_crs, dst_crs, xs, ys)
    return [[x,y] for x,y in zip(xs, ys)]

with rasterio.open('/tmp/b1_moll.tif') as dataset:
    src_crs = 'EPSG:4326'
    dst_crs = dataset.crs.to_proj4()  # '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs'
    # dst_crs = dataset.crs.to_wkt()  # 'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS    84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

coords = [[123.456, -34.567]]  # [longitude, latitude] not [lat, lon]...
    new_coords = reproject_coords(src_crs, dst_crs, coords)

values = list(rasterio.sample.sample_gen(dataset, new_coords))

for (lon, lat), value in zip(coords, values):
        print(lon, lat, value[0])  # value[0] == band 1 value at lon, lat

Ivotje50 · Answer

The above solution used to work pretty good, but the transformation of coordinates from one system to another has been modified (I think Fiona doesn't even include it anymore).
Best thing is to use pyproj Transformer, as shown below.
import rasterio as rio
from pyproj import Transformer
with rio.open(tif) as dataset:
    examples = [(52.2271, 4.9013), (52.3778, 4.8816)]
    transformer = Transformer.from_crs("epsg:4326", dataset.crs)
    coords = [transformer.transform(x, y) for x,y in examples]

Read more on how to best do this here: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1. Hope this helps!

Using Rasterio to index by traditional longitude and latitude with Molleweide Projection GeoTIFF

2 Answers

Add your own answers!

Ask a Question