Rasterio's features.rasterize()
function is not working as expected.
I have two files:
- OpenStreetMap shapefile with land usages
- 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
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.