3

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

2

2 Answers 2

3

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
2
  • This should be the solution. The dataset did not match my expectations however.
    – Ewin Zuo
    Mar 28, 2020 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, 2020 at 3:56
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!

3

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.