Building a Rasterio Workflow for Your Project

Raster Metadata

Code Block 1: dtype differences between line 9 vs line 13 are highlighted in red
  • Line 4: We load the raster into an array but force it to be read as a Float32. We do this by adding .astype(rasterio.float32) after f.read(). You can call all the datatypes such as rasterio.uint16, rasterio.int8, rasterio.float64, and more. If your IDE supports auto-completion you can see the other options.
  • Line 5: We load the metadata. The metadata is a Python dictionary object so we can manipulate as such.
  • Line 7: The raster is loaded into a standard NumPy array as seen with the print statement. So you can treat it as an ordinary NumPy array and even use NumPy functions with it.
  • Line 8–9: This is the important part. Even though we force Rasterio to read the raster as a Float32 the metadata does not follow it and it still says that the image is a uint16.
  • Line 11: We update the metadata dictionary using the dtype property of our Raster array.
  • Line 13: The metadata is now updated which is highlighted in red for clarity.
Code Block 2
  • Line 1: This is our metadata from code block 1. So we see we have our modified dtype value.
  • Line 2: This shows that the raster is also a Float32 so we verify that the dtypes in the array and metadata match so we can proceed to create our new raster.
  • Line 3: Similar to your standard Python context manager. The metadata has to be used with the double asterisks “**” which is a special syntax for Python arguments.

Sample Workflow: Calculating NDVI

Image from Digital Globe
Code Block 3
  • Line 5: We anticipate floating point calculations so we read the raster as a Float32.
  • Line 11–12: Band order is Blue, Green, Red, Near-Infrared. Thus, using Python zero-based indexing we extract the red band by using 2 and the NIR band using 3.
  • Line 14: This is the formula for NDVI.
  • Lines 15–17: The output of the NDVI calculation is a flat 2 dimensional array. Rasterio still needs a 3 dimensional array even if it is a single band image. Remember that this is still a NumPy object so we can use np.expand_dims to add an extra dimension to the zero axis.
  • Lines 19–20: NDVI is a single band image so we want to update the metadata to match this. This is under the count property. A good practice is to assign it to the first index of the shape property.
  • Line 22: We can see our updated metadata for our final NDVI output.

Final Thoughts

--

--

--

Researcher working on geospatial sciences and general programming

Love podcasts or audiobooks? Learn on the go with our new app.

Recommended from Medium

CS 373 Spring 2022: Kyle Kamka

Druid Cluster Setup

4. Adding DynamoDB to Serverless Microservice

Testing For Cross Browser Compatible Web App using Joomla

DLithe_BC_NFS_T_Task14_WebDesigning

Stobox Weekly Report 05 (01.02–07.02)

Object Oriented Programming

How to Write Testable Code

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store
Panji Brotoisworo

Panji Brotoisworo

Researcher working on geospatial sciences and general programming

More from Medium

The Easiest Way to Plot Topography

Guidance of installing and usage of PostGIS on MacOS

Bringing Pulse data to life in real-world applications

Cloud GeoServer in 20 minutes