Last time I wrote about how Cloud Optimized Geotiffs (COGs) are a great way to save you from having to download data and that they can be be easily streamed into QGIS (3.4 and later). If you would like to read that post it is available here.
This time though I want to look in a bit more detail at accessing COGs using Python AND at how you can download a small subset of imagery using Python. Actually, this is one of the most common questions I have been asked.
“I don’t want to download a whole scene, just a few bands, over a small area. How can I do this?“
One of the less mentioned benefits of COGs is that you do have more control – by this I mean that you are not just downloading a zip file of one tile. In this post I will show you a simple way of working with COGs in Python. There is a really great guide that I have used as inspiration here. Please do check it out.
Cloud Optimized Geotiff – clipping and downloading part of a Landsat 8 scene step by step
First off, let’s import all the libraries we will be using.
%matplotlib inline import rasterio import matplotlib.pyplot as plt import numpy as np from rasterio.plot import show import geopandas as gpd from rasterio.mask import mask
I have already found the landsat 8 tile I want to stream – it’s here: http://landsat-pds.s3.amazonaws.com/c1/L8/202/025/LC08_L1TP_202025_20190420_20190507_01_T1/ I am going to get a list of the bands I want and order them so it’s descending. That way I can create a 432 RGB image.
landsatband =  landsat_fp = 'http://landsat-pds.s3.amazonaws.com/c1/L8/202/025/LC08_L1TP_202025_20190420_20190507_01_T1/LC08_L1TP_202025_20190420_20190507_01_T1_' bands=('4.TIF', '3.TIF', '2.TIF') for i in range(1,11): ## 11 bands of Landsat right? band = landsat_fp+'B'+str(i)+'.TIF' if band.endswith((bands)): landsatband.append(band) landsatband.sort(reverse=True) ## because I want the order 4,3,2 (but change if you need to!) print(landsatband)
If you want to add more bands to your image stack then add them into the tuple, bands. Next, I want to load the vector data containing the area of interest. I am going to open a geojson file with GeoPandas to do this:
import geopandas as gpd gdf = gpd.read_file('boundry.geojson')
Next, I will write a little function that I will call for each band to do the clipping. It’s called read_imagery and takes an image path as the only argument.
def read_imagery(path): data = rasterio.open(path) masked, mask_transform = mask(dataset=data, shapes=gdf.geometry, crop=True) print(masked.shape) return data, masked[0,:,:], mask_transform
If you want a detailed explanation of this have a look at my previous blog post here:
Now I will call this function as I iterate over my list of image paths. I’ll append each array to a list and then finally stack them into a 3d array called rgb.
lsarray= for path in landsatband: data, masked, mask_transform = read_imagery(path) lsarray.append(masked) rgb = np.stack(lsarray)
I can plot one band to check that it has worked using this line:
And then, finally, write it out to a tif:
profile = data.meta print(profile) WIDTH = rgb.shape ## get the dimensions of the image we are writting out HEIGHT = rgb.shape profile.update(count = 3, transform=mask_transform, height = HEIGHT, width = WIDTH) print(profile) with rasterio.open('clipCOG_432.tif', 'w', **profile) as dst: dst.write(rgb)
Again, I covered this in the above blog post. Open this image in QGIS (or any Geotiff reader):
And there you go. A simple guide to using a COG in Python to download an image subset based on an area you have defined. This is a simple example; sometimes your area of interest will cover two tiles and you’ll need to do some merging / mosaicing prior to clipping and downloading.
All the code (and some extra bits) is available here to download
I hope this has been of some help.
I am a freelancer able to help you with your projects. I offer consultancy, training and writing. I’d be delighted to hear from you. Please check out the books I have written on QGIS 3.4
I have grouped all my previous blogs (technical stuff / tutorials / opinions / ideas) at http://gis.acgeospatial.co.uk.
Feel free to connect or follow me; I am always keen to talk about Earth Observation.