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)]]}

Does this answer your question? Extracting data from a raster — snowman2, Sep 23 '20 at 01:20

user2856 · Accepted Answer · 2020-03-29T04:02:01.100

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

This should be the solution. The dataset did not match my expectations however. — Ewin Zuo, Mar 28 '20 at 22:28
@EwinZuo What do you mean. I just downloaded the raster, picked some coords in US & Australia and the results were correct. Don't forget coordinate order in GIS software is x,y i.e longitude, latitude not lat, lon. — user2856, Mar 29 '20 at 03:56

Ivotje50 · Answer 2 · 2020-09-23T09:10:33.570

3

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!

edited Sep 23 '20 at 09:10

answered Sep 22 '20 at 17:40

Ivotje50

281
2
7

I would avoid using to_proj4 https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems. I recommend instead passing in the rasterio.CRS object directly as it is supported by pyproj https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#converting-from-rasterio-crs-crs-to-pyproj-crs-crs. – snowman2 Sep 23 '20 at 01:05
See: https://gis.stackexchange.com/a/358058/144357 – snowman2 Sep 23 '20 at 01:19
@snowman2 looks good, I've modified it to feed the crs-object directly. Sometimes I'm surprised how Python "just works"! – Ivotje50 Sep 23 '20 at 09:10

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

2 Answers2