[GIS] Rasterio does not execute the rasterize function

fionagdalpythonrasterio

Rasterio's features.rasterize() function is not working as expected.

I have two files:

  1. OpenStreetMap shapefile with land usages
  2. Raster (.tif) with satellite imagery

Goal: Create a raster file from the OSM shapefile where attribute type = 'parking'.
Requirement: The raster file should have the same shape and projection as the .tif satellite imagery file.

Here is what they look like loaded into QGIS. You can see the feature polygons in green over the satellite background
they definitely overlap!

Here are the projection and transform

  • using rasterio.open().meta.

Raster file

{'affine': Affine(1.0, 0.0, 543673.0, 0.0, -1.0, 4192469.0),
 'count': 3,
 'crs': CRS({'init': u'epsg:26910'}),
 'driver': u'GTiff',
 'dtype': 'uint8',
 'height': 21509,
 'nodata': None,
 'transform': (543673.0, 1.0, 0.0, 4192469.0, 0.0, -1.0),
 'width': 17245}
  • and using fiona.open().meta

Vector file

{'crs': {'init': u'epsg:4326'},
'crs_wkt': u'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]',
'driver': u'ESRI Shapefile',
'schema': {'geometry': 'Polygon',
'properties': OrderedDict([(u'id', 'int:11'),
    (u'osm_id', 'float:19'),
    (u'name', 'str:104'),
    (u'type', 'str:17'),
    (u'area', 'float:32.10'),
    (u'z_order', 'int:11')])
}}

Here's the code:

import rasterio
import numpy as np
from rasterio import features
import fiona

#Load vector file with fiona
vector_file = 'shapefiles/sf_osm/san-francisco_california_osm_landusages.shp'

shapefile = fiona.open(vector_file)
geom = [shapes['geometry'] for shapes in shapefile]
attrib = [shapes['properties'] for shapes in shapefile]

#filter to features where type = parking
    type_match_idx = [i for i, feature in enumerate(attrib) if feature['type'] == 'parking']
    feature_geom = [geom[i] for i in type_match_idx]
    feature_attrib = [attrib[i] for i in type_match_idx]

#get metadata for satellite imagery
rst = rasterio.open('satellite.tif')

Everything is working fine up until here. This is where it stops acting as expected. 

#rasterize using the shape and transform of the satellite image
image = features.rasterize(feature_geom, out_shape=rst.shape, transform=rst.transform)

print(np.max(image)) #results in 0

This produced a numpy array with every value equal to zero. I would expect some to be non-zero where features exist.

#saving image
with rasterio.open(
        'rasterized-results2.tif', 'w',
        driver='GTiff',
        transform = rst.transform,
        dtype=rasterio.uint8,
        count=1,
        width=rst.width,
        height=rst.height) as dst:
    dst.write(image, indexes=1)

This produces a blank (all black) image.

Why is features.rasterize() producing no data?

Best Answer

The projections are different. They overlay in QGIS because QGIS reprojects on the fly.

Use fiona.transform.transform_geom to reproject the vector geometry from epsg:4326 to epsg:26910 before passing to rasterio.

Related Question