Digitizing in Python

Purpose

The purpose of this tutorial is to automate the digitization of land cover features such as vegetation, waterbodies, and roads, using Sentinel-2 satellite imagery and Python-based geospatial tools.

Objective

Use Google Earth Engine (GEE) to acquire the satellite images for processing in python. Compute NDVI (vegetation), NDWI (water), NDBI (roads) indices. Export rasters to your Google Drive. Convert rasters to vectors. Visualize and colourize digitized features.

Software Used

Google Earth Engine Google Colab or any other Python language system or notebooks.

Python libraries:

 earthengine-api, geemap, rasterio, geopandas, shapely, matplotlib

Setting up Colab and Google Earth Engine

Set up a Google Earth Engine account for free with these similar instructions (note that you do not need to make an account under a student requirement or the sort or any paid service, this is free!): https://courses.spatialthoughts.com/gee-sign-up.html

After your account is created, create a new project under a name of your choice. This will be linked back in the Colab notebook.

To set up a Colab notebook, simply open a new Colab notebook within your Google Drive.

Now, in the Colab notebook, follow this code and replace the project name with your own project name in the code snippet below:

 
!pip install earthengine-api geemap rasterio geopandas shapely matplotlib

import ee
import geemap
import geopandas as gpd
from shapely.geometry import shape
import matplotlib.pyplot as plt
import numpy as np
from rasterio.features import shapes
import rasterio
from rasterio.plot import show

ee.Authenticate()
ee.Initialize(project='your-project-name')

What is Digitizing?

Digitizing is the process of converting geographic features into digital vector forms. Though this used to be a manual process using different techniques by hand, it can be done in a digital space or even automated using indices and raster-to-vector conversion.

Why is Digitizing Important?

Digitizing geographic data is crucial for converting analog maps and spatial information into a digital space because digital formats are preferred, with the ability to be easily analyzed, shared and constantly up to date.

Data

• Sentinel-2 imagery from Google Earth Engine
• Area of Interest: Hogsback River near Carleton University

Choosing Your Own Satellite Image

The code below is based on the coordinate data of Hogsback River, near Carleton University. Any area in the world within the limits of what satellite images Google Earth Engine has access to can be used in place of the chosen coordinates for this tutorial.

 
aoi = ee.Geometry.Point([-75.6909, 45.3795])  # longitude, latitude
aoi_buffer = aoi.buffer(1000).bounds()        # 1 km buffer around the river

Methods

The main method is using NDVI, NDWI, and NDBI and setting up these bands to represent them in the output images as well as setting up a threshold value for each.

Spectral Indices

Index	Purpose	Bands Used	Threshold
NDVI	Vegetation	B8, B4	0.3
NDWI	Waterbodies	B3, B8	0.2
NDBI	Built-up Roads	B11, B8	0.2

Steps

Part 1

Once you have loaded the area of interest, the next step is to select the parameters for time, cloudiness percentage, and to see what the final output will look like.

Selecting filters for the image and library from Copernicus Sentinel-2 imagery:

 
collection = (ee.ImageCollection("COPERNICUS/S2_SR")
              .filterBounds(aoi_buffer)
              .filterDate("2024-06-01", "2024-08-31")  # summer imagery
              .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10)))

Select least cloudy image:

 
image = collection.sort("CLOUDY_PIXEL_PERCENTAGE").first()

Compute NDVI (this is just to see a third of what the output will resemble):

 
ndvi = image.normalizedDifference(["B8", "B4"]).rename("NDVI")

Selecting parameters to display the map:

 
Map = geemap.Map()
Map.centerObject(aoi_buffer, 13)
Map.addLayer(image, {"bands":["B4","B3","B2"], "min":0, "max":3000}, "Sentinel-2 RGB")
Map.addLayer(ndvi, {"min":0, "max":1, "palette":["white","green"]}, "NDVI")

Saving the raster to the Google Drive:

 
task = ee.batch.Export.image.toDrive(
    image=ndvi,
    description="Hogsback_NDVI",
    folder="EarthEngine",
    fileNamePrefix="hogsback_ndvi",
    region=aoi_buffer.getInfo()['coordinates'],
    scale=10,
    crs="EPSG:4326"
)
task.start()
print("Export started. Check Google Drive → EarthEngine folder for hogsback_ndvi.tif")

Mount Google Drive in Colab once the export finishes:

 
from google.colab import drive
drive.mount('/content/drive')

Load the raster locally with Rasterio. (Note that the filepath will look different depending on where you have it saved in your drive):

dataset = rasterio.open("/content/drive/My Drive/EarthEngine/hogsback_ndvi.tif")
ndvi_array = dataset.read(1)