Apply same coordinate system to raster image and GeoJSON with rasterio

Question

I have a dataset containing raster images and corresponding features in geojson files. I would like to use rasterio.features.rasterize to create a footprint image of the features on the raster image, but the resulting footprint array stays empty. I think the issue is that the raster image and the features are defined for different coordinate systems.

The raster image uses EPSG-3857:

>>> data = rio.open('image.tif')
>>> data.crs
CRS.from_epsg(3857)

and the underlying transformation seems to be defined in bounding box coordinates:

>>> data.transform
Affine(4.77731426716, 0.0, -1942112.014669758,
       0.0, -4.77731426716, 1658377.7656751834)

The geojson file on the other hand seems to be defined using WGS 84 coordinates:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "Id": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -17.441005042849905, 14.689811850766825 ], [ -17.441172607380693, 14.689811850766825 ], [ -17.441172607380693, 14.689967040094865 ], [ -17.441005042849905, 14.689967040094865 ], [ -17.441005042849905, 14.689811850766825 ] ] ] } },
...

If I apply the affine transformation from the image to pixel coordinates, I end up with coordinates in bounding box coordinates:

>>> data.bounds
BoundingBox(left=-1942112.014669758, bottom=1653485.7958656116, right=-1937220.0448601863, top=1658377.7656751834)
>>> data.transform * (100, 100)
(-1941634.283243042, 1657900.0342484673)

My two questions are:

1. How do I transform these coordinates into world coordinates? And...
2. Is there a way to bring the geojson file coordinates on the same system as the raster image to generate the footprint?

As a reference, I wrote the following script to rasterize the features on the image:

>>> import rasterio as rio
>>> from rasterio import features
>>> import geopandas

>>> data = rio.open('image.tif')
>>> df = geopandas.read_file('labels.geojson')

>>> feature_list = list(zip(df['geometry'], [255]*len(df)))
>>> fpt = features.rasterize(shapes=feature_list,
                             out_shape=data.read(1).shape,
                             transform=data.transform)

>>> print(fpt.sum()) # to show that the footprint raster is empty
0

Answer 1

As a reference for myself, or anybody else who might have this problem in the future:

(1.) reproject raster data to crs used by features (EPSG:4326 in this case)

import numpy as np
import rasterio as rio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio import features
import geopandas

dst_crs = 'EPSG:4326'  # WGS 84

df = geopandas.read_file('labels.geojson')
feature_list = list(zip(df['geometry'], [255]*len(df)))

with rio.open('images.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    srcdata = np.array([src.read(i) for i in src.indexes])
    dstdata = np.zeros(srcdata.shape)
    
    reproject(source=srcdata,
              destination=dstdata,
              src_transform=src.transform,
              src_crs=src.crs,
              dst_transform=transform,
              dst_crs=dst_crs,
              resampling=Resampling.nearest)

coo = transform * (100, 10) # coordinate transformation to WGS 84  

(2.) transform features to common crs:

import rasterio as rio
from rasterio import features
import geopandas

data = rio.open('image.tif')
df = geopandas.read_file('labels.geojson')

allfeatures = [f['geometry'] for f in df.__geo_interface__['features']]
allfeatures_reprojected = rio.warp.transform_geom("EPSG:4326",
            data.crs, allfeatures)

feature_list = list(zip(allfeatures_reprojected,
                        [255]*len(allfeatures)))

fpt = features.rasterize(shapes=feature_list,
                         out_shape=data.read(1).shape,
                         transform=data.transform)

Please note that the latter example requires rasterio>=1.2dev.