1

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
1
  • 2
    Welcome to GIS.SE. Please focus only on one question, otherwise answering and later on for other users finding the correct/relevant answer is unnecessarily hard.
    – Erik
    Aug 11 '20 at 8:06
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.