How to Plot Coordinates on Landsat Satellite Images with Python | by Conor O’Sullivan | May, 2023

By Jessie Hobb On May 17, 2023

Mapping pixel positions to geographic coordinates using Landsat metadata and Rasterio

Location, location, location! Not just a housing market cliché but also incredibly important to remote sensing. Whether it’s for monitoring environmental changes, analyzing urban development, or tracking crop health, precise geolocation is crucial. We need to know the exact coordinates of objects in satellite images to ensure accurate analysis and interpretation.

So, we’ll explore how to plot coordinates directly onto Landsat scenes using two approaches:

The Landsat metadata file (MLT)
Rasterio — a package used to access geospatial raster data

We’ll also use Rasterio to reproject the coordinates of a satellite image. Specially, we’ll go from the original coordinate reference system (UTM) to the one used by google maps (EPSG:4326). Along the way, we’ll discuss the code and you can find the full project on GitHub.

We need to start by downloading a Landsat scene. You can do this using the EarthExplorer portal. Alternatively, if you want to use Python, the article below takes you through the process:

In the end, you should have a folder similar to Figure 1. These are all the files available for a Landsat level 2 science product. We’ll be working with the red visible band (*_B4.TIF) and the JSON metadata file (*_MTL.json). These are highlighted below.

A list of Landsat level-2 science product files. The red visible light band and MLT metadata file have been highlighted. — Figure 1: Landsat level-2 science product files (source: author)

This particular scene was taken above the east coast of Ireland. To see this we visualise the red band. We load the tiff file (line 9). We scale the band (line 12) and clip it to enhance it’s contrast (line 15). You can see the resulting image in Figure 2.

import tifffile as tiff
import numpy as np
import matplotlib.pyplot as plt# Landsat scene ID 
ID = 'LC08_L2SP_206023_20221118_20221128_02_T1'
# Load Red (B4) band
B4 = tiff.imread('./data/{}/{}_SR_B4.TIF'.format(ID, ID))
# Scale band
B4 = np.clip(B4*0.0000275-0.2, 0, 1)
# Clip to enhance contrast
B4 = np.clip(B4,0,0.2)/0.2
# Display grayscale image
fig, ax = plt.subplots(figsize=(10, 10))
plt.imshow(B4, cmap='gray')
ax.set_axis_off()

Visualisation of the red visible light band of Landsat scene of the east coast of Ireland — Figure 2: Visible red light band of Landsat scene of the east coast of Ireland (source: author)

We want to plot a point for Howth, an Irish village, directly onto this image. Looking at Figure 2, you can see it sticking out at about 1/3 of the way down the coast. The specific latitude and longitude points we want to plot are given in Figure 3.

Google Maps screenshot of the latitude and longitude coordinates for Howth, Dublin — Figure 3: Latitude and longitude coordinates for Howth, Dublin (source: google maps)

Okay, so plotting points with the MLT file is not as straightforward as using Rasterio. Still, it is worth exploring this method. You will see that we do not have to rely on a package. We will also learn about the coordinate reference system (CRS) used by Landsat.

Landsat metadata

The MLT file contains a lot of metadata about the Landsat scene. You can find all the details in the data format control book. We are interested in the projection attributes. These are details about the scenes CRS and will allow us to geolocate the pixels in Figure 2. We load these attributes (lines 4–5).

import json# Load metadata file
MTL = json.load(open('./data/{}/{}_MTL.json'.format(ID, ID)))
projection = MTL['LANDSAT_METADATA_FILE']['PROJECTION_ATTRIBUTES']
projection

Figure 4, gives the list of projection attributes. The first item tells us that UTM is the CRS for this scene. This is a projected coordinate system (PCS). It is used to project geolocation on the Earth surface (a sphere) onto the scene (2D surface). Later we will see how this differs from a geographic coordinate system (GCS).

A list of projection attribute for a Landsat scene taken from the MLT file — Figure 4: projection attribute for Landsat scenes (source: author)

When converting UTM to pixel positions, we will need two other pieces of information from Figure 4. These are:

GRID_CELL_SIZE_REFLECTVE — this tells us that each pixel covers 30 square meters of land
Projection coordinates (e.g. CORNER_UL_PROJECTION_X_PRODUCT) — these give the UTM coordinates (in meters) for the scenes bounding box

In fact, we only need the first 2 coordinates. These are for the upper left (UL) corner of the bounding box. Looking at Figure 5, we see how this information is used to convert to pixel position. We know each pixel is 30m². This means that we can start at (543000, 6004200) and add/subtract 30m to get the start of the next pixel. We select a pixel if the UTM coordinates fall in that pixel’s range.

The UTM coordinates for start and end of pixels taken from the upper left corner of a landsat image. — Figure 5: UTM coordinates of start and end of pixels

Plotting coordinates on a Landsat band

We use this information to create the coords_to_pixels function. As parameters, it takes the UTM coordinates for the UL corner (ul) and the coordinate we want to convert (utm).

def coords_to_pixels(ul, utm, m=30):
""" Convert UTM coordinates to pixel coordinates"""x = int((utm[0] - ul[0])/m)
y = int((ul[1] - utm[1])/m)
return x, y

We can now use this function to find the pixels for our point on Howth. We start by getting the UTM coordinates for the UL corner (lines 4–5). We took our points straight from google maps and so they are given in latitude and longitude (lines 8–9). We convert these to UTM coordinates (line 12) which gives us values of (694624, 5918089). Finally, we convert these to pixels (line 16). This gives us pixel values of (5054, 2870).

import utm # Get UTM coordinates of upper left corner
ul = [float(projection['CORNER_UL_PROJECTION_X_PRODUCT']),
float(projection['CORNER_UL_PROJECTION_Y_PRODUCT'])]
# Lat/long coordinates of Howth
lat=53.3760
long=-6.0741
# Convert GPS to UTM coordinates
utmx,utmy,_,_ = utm.from_latlon(lat,long)
print(utmx,utmy)
# Convert UTM to pixel coordinates
x,y = coords_to_pixels(ul, [utmx,utmy])
print(x,y)

Let’s see how accurate this is. We add a white circle to the image in Figure 2 at the given pixel values (line 5). To see our point we zoom in by cropping the image around the circle (lines 8–9). Comparing Figure 6 to the point in Figure 3, we can see that we have accurately plotted the point onto the Landsat scene. In general, 90% of all Landsat points will have a positional error of less than 12 meters.

import cv2 # Add circle to image at pixel coordinates
image = B4.copy()
image = cv2.circle(image, (int(x), int(y)), 10, 1,5)
# Crop image around pixel coordinates
h = 100
crop_image = image[y-h:y+h, x-h:x+h]

White circle plotted on a lansat image which give the coordinates for Howth — Figure 6: Howth coordinates plotted using Landsat Metadata file (source: author)

The above process works! But, we can make our lives easier using the Rasterio package. We open the tiff file for the red visible light band (line 4). Metadata about the scene is embedded in the file. For example, we print the CRS (line 7) which outputs “EPSG:32629”. This is the EPSG code for UTM zone 29N. In other words, the same CRS that was given in the MLT file.

import rasterio as rio# Open red channel with rasterio
B4 = rio.open('./data/{}/{}_SR_B4.TIF'.format(ID, ID))
# Display image map projection 
print(B4.crs)

We can use the index function to convert UTM to pixel coordinates (line 9). This outputs (5054, 2870). The same as before!

# Lat/Long coordinates of Howth
lat=53.3760
long=-6.0741# Get UTM coordinates
utmx,utmy,_,_ = utm.from_latlon(lat,long)
# Convert UTM to pixel coordinates
y,x = B4.index(utmx,utmy)
print(x,y)

So, with Rasterio we no longer have to worry about the bounding box or even the CRS. However, we still have to convert latitude and longitude to UTM. If we want to plot the latitude and longitude directly, we will have to reproject the Landsat scene to the Google Maps CRS.

Reprojecting a Landsat scene

We do this using the code below. In summary, it takes the red band image, reprojects it to EPSG:4326 (Google Maps CRS), and saves the reprojected image to a new file. Some things to point out:

We calculate a new transform function (line 12). This is what Rasterio uses to convert coordinates to pixels.
When we reproject the file we use the nearest resampling method (line 42). This means we use the value of the nearest cells to fill in the values of other cells.

We need a resampling method because the reprojection process will “warp” the image to fit the new CRS. We will see this when we output the reprojected image.

from rasterio.warp import reproject, Resampling, calculate_default_transformdst_crs = "EPSG:4326"  # google maps CRS
filename = './data/{}/{}_SR_B4.TIF'.format(ID, ID)
new_filename = './data/{}/{}_SR_B4_EPSG4326.tif'.format(ID, ID)
# Open file
with rio.open(filename) as src:
src_transform = src.transform
# calculate the transform matrix for the output
dst_transform, width, height = calculate_default_transform(
src.crs,
dst_crs,
src.width,
src.height,
*src.bounds,  # unpacks outer boundaries
)
# set properties for output
dst_kwargs = src.meta.copy()
dst_kwargs.update(
{
"crs": dst_crs,
"transform": dst_transform,
"width": width,
"height": height,
"nodata": 0,  
}
)
# write to disk
with rio.open(new_filename, "w", **dst_kwargs) as dst:
# reproject to new CRS
reproject(
source=rio.band(src, 1),
destination=rio.band(dst, 1),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)

To confirm the process has worked, we load the projected image (line 2) and output its CRS (line 3). This gives us “EPSG:4326”. Success!

# Open red channel
B4 = rio.open('./data/{}/{}_SR_B4_EPSG4326.TIF'.format(ID, ID))
print(B4.crs)

We can now get the pixel coordinates directly from lat/long coordinates. This gives a value of (6336, 2227). These are different to before because, as mentioned, the projected image has been warped.

# GPS coordinates of Howth
lat=53.3760
long=-6.0741# Convert UTM to pixel coordinates
y,x = B4.index(long,lat)
print(x,y)

We can see this when we plot the point on the projected image in Figure 7. Comparing it to Figure 6, the aspect ratio has changed. In fact, many of the pixels will have been altered. This is because, unlike UTM, EPSG:4326 is a geographic coordinate system (GCS). It is used to plot points on a sphere (i.e. the Earth). Trying to display an image with this CRS on a 2D plane leads to the distorted aspect ratio.

White circle plotted on a lansat image which give the coordinates for Howth. The image is warped due to the reprojection of the scene’s CRS. — Figure 6: Howth coordinates plotted on the projected image (source: author)

This article has focused on plotting coordinates on a single band. If you are working with multiple bands or RGB visualisations, the underlying method will be the same. This is because all bands in the Landsat scene will have the same CRS. See the article below for details on visualising the RGB channels.