Forem: Kofi Adu Agyekum

Calculating NDVI in Python with Rasterio and GeoPandas

Kofi Adu Agyekum — Wed, 04 Jun 2025 07:11:36 +0000

After learning a lot from my first attempt at NDVI processing (where just about everything went wrong), I built a clean, working NDVI pipeline using Python. This post walks through how I calculated and visualized NDVI over the Pohjois-Savo region in Finland using Landsat 8 data.

What is NDVI?

NDVI stands for Normalized Difference Vegetation Index. It’s a widely used metric in remote sensing that helps identify healthy vegetation using the difference between red and near-infrared (NIR) light reflected by surfaces.

The formula:

NDVI = (NIR - Red) / (NIR + Red)

Higher NDVI values generally indicate denser, healthier vegetation.

Tools & Data

Landsat 8 Surface Reflectance data from USGS EarthExplorer
- Band 4 = Red
- Band 5 = Near Infrared (NIR)
Python Libraries:
- rasterio (for raster data)
- geopandas (for shapefiles/geojson)
- numpy (for array math)
- matplotlib (for visualization)
Boundary: A GeoJSON file covering part of the Pohjois-Savo region in Finland

Step-by-Step Workflow

1. Import libraries

import rasterio
import geopandas as gpd
import numpy as  np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from rasterio.mask import mask
from rasterio.crs import CRS
from rasterio.features import shapes
from shapely.geometry import shape

2. Load the Red, NIR Bands, and Reproject the GeoJSON Boundary

red_path = '../data/raw/band4.TIF'
nir_path = '../data/raw/band5.TIF'
shape_path = '../data/raw/workarea.geojson'

3. Load, Clip the Bands to the Boundary and Calculate NDVI

#open bands
with rasterio.open(red_path) as red_src, rasterio.open(nir_path) as nir_src:
    red = red_src.read(1).astype('float32')
    nir = nir_src.read(1).astype('float32')

    #check for invalid crs (EPSG:9001 or None)
    def fix_crs(src):
        if not src.crs or "9001" in str(src.crs):
            return CRS.from_epsg(32635)  
        return src.crs

    red_crs = fix_crs(red_src)
    nir_crs = fix_crs(nir_src)

    #use red's crs as base
    working_crs = red_crs

    #load and reproject farm boundary
    #if None you can set to 3857 also   
    farm_shape = gpd.read_file(shape_path)
    if farm_shape.crs is None:
        farm_shape.set_crs("EPSG:4326", inplace=True)
    farm_shape = farm_shape.to_crs(working_crs)

    #clip raster bands
    clipped_red, _ = mask(red_src, farm_shape.geometry, crop=True)
    clipped_nir, _ = mask(nir_src, farm_shape.geometry, crop=True)

    #recalculate NDVI
    ndvi_clipped = (clipped_nir[0] - clipped_red[0]) / (clipped_nir[0] + clipped_red[0])

4. Plot your graph

plt.figure(figsize=(8, 6))
plt.imshow((nir - red) / (nir + red), cmap='RdYlGn', vmin=-1, vmax=1)
plt.colorbar()
plt.title("Full NDVI")
plt.axis("off")
plt.show()

You will see something like this:

5. Classify NDVI (Optional)

def classify_ndvi(ndvi_array):
    classified = np.full(ndvi_array.shape, -1, dtype='int8') #-1 = no data
    classified[(ndvi_array >= -1) & (ndvi_array < 0.2)] = 0  #barren/water
    classified[(ndvi_array >= 0.2) & (ndvi_array <= 0.5)] = 1 #moderate vegetation
    classified[(ndvi_array > 0.5) & (ndvi_array <= 1)] = 2 #healthy/high vegetation
    return classified

# appling to ndvi array
classified_ndvi = classify_ndvi(ndvi_clipped)

6. Plot Classified Graph

colors = ['skyblue', 'yellowgreen', 'darkgreen']
cmap = mcolors.ListedColormap(colors)
labels = ['Barren/Water', 'Moderate', 'Healthy']

plt.figure(figsize=(8, 6))
img = plt.imshow(classified_ndvi, cmap=cmap, vmin=0, vmax=3)
cbar = plt.colorbar(img, ticks=[0.5, 1.5, 2.5])
cbar.ax.set_yticklabels(labels)
plt.title("NDVI Vegetation Classification")
plt.axis("off")
plt.tight_layout()
plt.show()

The result will look something like this:

7. Polygonise your data

with rasterio.open(red_path) as src:
    transform = src.transform
    crs = src.crs

#generate dictionary (geometry, value) from classified ndvi
polygons = []
for geom, value in shapes(classified_ndvi, transform=transform):
    if value != -1:#skips no data
        polygons.append({"geometry": shape(geom), "value": int(value)})


#convert to geodataframe
gdf = gpd.GeoDataFrame(polygons, crs=crs)

#adding labels
gdf["label"] = gdf["value"].map({0: "Barren", 1: "Moderate", 2: "Healthy"})

#save as geojson
gdf.to_file("classified_ndvi_polygons.geojson", driver="GeoJSON")

#preview
gdf.head()

8. GeoJson, Shapefiles and GeoPackage

#save as geojson file
gdf.to_file("classified_ndvi_polygons.geojson", driver="GeoJSON")

#save as a .shp file
gdf.to_file("classified_ndvi_polygons.shp")

#save as a geopackage
gdf.to_file("classified_ndvi.gpkg", driver="GPKG")

While working on this project, I also ran into a few format-related issues when exporting the classified results. At first, I saved the vectorized NDVI zones as a GeoJSON, but the file ended up being around 352MB, which was too large for QGIS to handle smoothly. It either failed to load or froze entirely. I tried again using a Shapefile, which worked better for loading large datasets, but came with limitations like short field names and no support for special characters. Finally, I exported the same data as a GeoPackage (.gpkg), and it worked perfectly. It handled the large file size, preserved full attribute names, and opened cleanly in QGIS. Based on this, I’d recommend GeoJSON for lightweight data, Shapefiles for compatibility, and GeoPackage for larger or more complex outputs.

9. Loading the Result in QGIS

Once I had the classified NDVI zones vectorized and saved as a GeoPackage, I brought the file into QGIS to generate the final map layout.

QGIS handled the GeoPackage smoothly — no freezing or lag, even with thousands of polygons. I applied a simple categorized style based on NDVI class (barren, moderate, healthy) for easy interpretation.

From there, I added:

A base map for reference (using XYZ Tiles)
A legend, scale bar, and labels
A title and layout design using the Print Layout tool

Then I exported the result as a high-resolution PDF map.

This step wasn’t just about visualization — it was about validating that the Python-generated output aligned with the real-world landscape and could be presented clearly.

Result

The final NDVI map showed strong vegetation density in forested areas and lower NDVI in more open or cleared zones. I processed everything in Python and confirmed that the output matched what I saw later in QGIS when visualizing the same area.

Takeaways

Always check CRS and bit depth before processing
NDVI outputs may need classification to be more interpretable
Vectorizing large rasters can produce heavy files — use GeoPackage or Shapefile for big outputs

If you want to see the mistakes that led to this working version, check out:

Before I Got the Map, I Got It All Wrong

Before I Got the Map, I Got It All Wrong

Kofi Adu Agyekum — Wed, 21 May 2025 10:59:40 +0000

I set out to calculate NDVI and visualize vegetation in Finland using satellite imagery. But before I could get anything meaningful on the screen, everything broke, and every fix taught me something new.

This post is a breakdown of what I got wrong, what didn’t make sense, and how I eventually got it working.

1. I Downloaded the Wrong Data Over and Over

I started by downloading Sentinel-2 .tiff files from Copernicus. They looked fine, no errors, so I loaded them into QGIS.

What I saw: a completely dark image.

I assumed the darkness was normal because it was “raw data.” But something didn’t feel right.

I deleted the files and downloaded another set, same issue. I tried different pixels, formats, even downloaded 32-bit, then 16-bit versions.

Still dark.

2. CRS Confusion Made It Worse

When I first loaded the raster in QGIS, the project was set to EPSG:3857 (Web Mercator). The raster wasn’t distorted — it just appeared much larger in scale compared to the base map. It didn’t fit well within the map canvas, and the proportions felt off.

Thinking the issue was related to projection, I changed the project CRS to EPSG:4326. But that actually made things worse — the raster became distorted, and the base map was visibly stretched.

That’s when I realized something important:

My raster’s pixel values were too high, and it couldn’t display correctly over the default basemap.

3. 8-Bit Fixed the Fit, Finally!

I went back and downloaded the 8-bit version of the data — and suddenly, it snapped into place on the map.

For the first time, it aligned properly with the basemap in QGIS.

That’s when it clicked: the issue wasn’t just the CRS, it was the raster’s bit depth and how QGIS handles rendering.

High-bit-depth rasters (16-bit or 32-bit) are designed for analysis, not visualization. They carry more numerical precision but often appear black, washed out, or oversized when layered over a map — especially without stretching or resampling.

Lesson learned:

If you’re doing visual mapping or quick previews, use 8-bit data.

If you’re doing NDVI, classification, or scientific analysis — 16-bit or 32-bit is fine, and you don’t even need a basemap.

Not everything needs to look good on a map — especially when the goal is analysis.

4. NDVI Calculation Gave Me... a Black Block

Now confident I had good data, I calculated NDVI in Python:

ndvi = (nir - red) / (nir + red)

Plotted it, and saw nothing but a solid black rectangle.

No zones, no variation, just... black.

5. Identify Tool Showed Nothing

I opened the NDVI result in QGIS and used the Identify Tool.

The values? Completely empty.

I ran the same raster through a Python check:

print(np.nanmin(ndvi), np.nanmax(ndvi))

→ Both returned NaN.

The entire image was full of NaNs.

At that point, I had tried almost every combination of download and reprocessing. Nothing was working.

6. Clipping and Vectorizing Worked, But the Result Was Weird

Next, I clipped the NDVI raster using a farm boundary and then vectorized the classified result (low, medium, high vegetation).

It worked — technically — but the output looked strange.

In QGIS, I saw massive red, yellow, and green blocks. The low NDVI zone appeared as a huge red shape across the map.

It didn’t look wrong, but it didn’t look quite right either.

7. I Asked for Help And Got the Real Data

Finally, I posted a quick question on r/gis:

The responses pointed me to USGS EarthExplorer, and this changed everything.

I downloaded Landsat 8 Level-2 Surface Reflectance, Bands 4 and 5 — full-resolution, proper metadata, 80–100MB each.

I re-ran my workflow. Everything clicked:

QGIS rendered the bands clearly
Python returned actual NDVI values
The visualization was clean and detailed

8. Saving as GeoJSON Was the Final Trap

With a great result in hand, I decided to export the classified NDVI zones to GeoJSON.

gdf.to_file("classified_ndvi.geojson", driver="GeoJSON")

It saved, but the file was 352MB.

When I tried to load it into QGIS, it froze. Or loaded partially. Or refused to render at all.

9. Shapefiles and GeoPackages Saved the Day

After some quick Googling, I realized GeoJSON doesn’t handle large vector data well.

So I re-exported the same data:

# Save as Shapefile
gdf.to_file("classified_ndvi.shp")

# Save as GeoPackage
gdf.to_file("classified_ndvi.gpkg", driver="GPKG")

Both formats worked flawlessly. QGIS loaded them in seconds, with smooth rendering and no crashing.

Final Thought

This was my first real geospatial project. I thought it would be simple. It wasn’t but that’s what made it valuable.

Every “wrong” step helped me understand something no tutorial could fully explain.

If you're working with satellite data: check your CRS, validate your inputs, and don’t be afraid to start over.

And when in doubt? Ask a community. A single Reddit reply saved this entire project.

Getting started with Openlayers in React

Kofi Adu Agyekum — Fri, 29 Mar 2024 19:43:39 +0000

CONTENT

Introduction to Maps and Openlayers
Openlayers
Attributes
Setting up project
Summary

Introduction to Maps and OpenLayers

Ever dreamt of building interactive maps to explore or showcase data? This article dives into the powerful duo of React and OpenLayers, transforming your web pages into dynamic geo-discovery tools.

Openlayers

Openlayers is a JavaScript library that is used to render onto or embed maps into a webpage. It is a free and opensource library that allows users to display maps and load vector data and markers from any data source. (https://openlayers.org/)

It is a powerful and versatile tool that allows developers to build a range of applications, from simple maps to complex, and highly interactive geospatial applications with JavaScript.

Attributes

Openlayers has 4 key components that are used to create and render maps. These are Map, View, Source and Layer.

Map, in openlayers, is the map container for rendering any map. Imported from the module: ol/map , it needs a target container (e.g a div) to contain the map.



import Map from "ol/map"

const map = new Map({target: 'map'})

To display the map on your webpage, simply add an 'id' attribute to the container element of the map with the same name as the target in the newly created map above. Feel free to choose any name for the target as long as it matches the HTML element's ID.



<div id="map" style="width: 100%; height: 400px"></div>

View is the part of the map creation responsible for centering, zooming and projection of the map. It is imported from the ol/view module. A View comes with a projection, which decides how coordinates are handled and what units are used for map resolutions. If you don't choose one, it defaults to Spherical Mercator (EPSG:3857), using meters as the unit for maps. The view attribute is added as an option to the new map created in the Map column above. ```tsx

import Map from "ol/map"
import View from "ol/view"

const map = new Map({
target: 'map',
view: new View({
center: [0, 0],
zoom: 2,
}),
})

The center attribute specifies the initial center coordinate of the map view. It's an array with two elements representing the longitude and latitude of the center point on the map. In the above example, the map will be centered where both the longitude and latitude is zero (0). 

The zoom attribute sets the initial zoom level of the map. It is used to set how magnified or shrunk the map is compared to its default setting. You also have the option to define both _minZoom_ and _maxZoom_, allowing you to specify the minimum and maximum levels of zoom that users can apply to the map.

For example: let's make the Colosseum in Rome our starting point. Just plug in the EPSG:3857 coordinates, which you can find here: https://epsg.io/map#srs=4326&x=0.000000&y=0.000000&z=1&layer=streets. Then, zoom in or out as you please!

```tsx


import Map from "ol/map"
import View from "ol/view"

const map = new Map({
  target: 'map', 
  view: new View({
      center: [1390659.798668, 5144570.023792],
      zoom: 17,
  }),
})

The result will look like this:

As you can see, the map is centered on the Colosseum.

Source is used to define the origin of spatial data that is displayed on the map. Openlayers support a variety of sources such as Tile sources like OpenStreetMap (OSM), Image sources like Web Map Service (WMS), Vector sources and XYZ sources.



import OSM from 'ol/source/osm';

const source = OSM();

This code grabs the OpenStreetMap data and prepares it to be shown on a map.

Layer, simply put is the visual representation of a data source. Examples include Tile layers, Image Layers and Vector Layers.



import TileLayer from 'ol/layer/tile';

const layer = new TileLayer({source: source});

For further reading visit (https://openlayers.org/doc/tutorials/concepts.html)

Setting up Openlayers in React

To create a map, you can merge all the essential components into one script as shown below:

First set up openlayers in react:



npm install ol
# or
yarn add ol

In the MapView component (Typescript):



import React, { useEffect } from "react";
import "ol/ol.css";
import Map from "ol/map";
import View from "ol/view";
import TileLayer from "ol/layer/tile";
import OSM from "ol/source/osm";

function MapView() {
  useEffect(() => {
    const map = new Map({
      target: "map",
      layers: [
        new TileLayer({
          source: new OSM(),
        }),
      ],
      view: new View({
        center: [0, 0],
        zoom: 2,
      }),
    });

    return () => {
      map.setTarget(null);
    };
  }, []);
return <div id="map" style={{width: "100%", height: "400px"}}/>;
}

export default MapView

The return statement in the useEffect ensures proper cleanup and management of resources associated with the map, thereby maintaining system efficiency and stability.

MapView component (Javascript):



import React, { useEffect } from "react";
import Map from "ol/Map.js";
import View from "ol/View.js";
import TileLayer from "ol/layer/Tile.js";
import OSM from "ol/source/OSM";
import "ol/ol.css";

function MapView() {
  useEffect(() => {
    const map = new Map({
      target: "map",
      layers: [
        new TileLayer({
          source: new OSM(),
        }),
      ],
      view: new View({
        center: [0, 0],
        zoom: 2,
      }),
    });

    return () => {
      map.setTarget(null);
    };
  }, []);

  return <div id="map" style={{ width: "100%", height: "400px" }} />;
}

export default MapView;

The above code snippet will display a basic OpenStreetMap (OSM) like this:

Summary

In our review, we've explored OpenLayers, a tool designed for creating maps, and delved into its integration with React. Through this exploration, we've gained insights into showcasing maps on our web pages by combining OpenLayers with React.