Rasterio GeoTIFF Indexing – Using Traditional Longitude and Latitude with Molleweide Projection

geotiff-tiffpythonrasterrasterio

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

Best 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
Related Question