Introduction
In this tutorial, you will do the following:
- Mosaic Basemap Quads
- Clip to an AOI
- Perform NDVI on the clipped mosaic
Before Steps Before we dive into this tutorial, you will need to have downloaded at least 2 quads from an area of your choice. In order to do this, please refer to this tutorial here for using the GUI or here for using API.
You will also need to have a polygon feature class for your area of interest which we will use for the clipping section. If you do not have ready, please use the following website to create one easily.
Along with the above, you will also need GDAL (Geospatial Data Abstraction Library) and its python bindings installed to run the commands below
1. Mosaic Basemap Quads¶
After downloading the quads from an area of your choice, running the following command will give you the spatial metadata of the quads.
!gdalinfo data/L15-0697E-0963N.tif
The response should look like below:
Driver: GTiff/GeoTIFF
Files: L15-0697E-0963N.tif
Size is 4096, 4096
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
Origin = (-6398696.510913709178567,-1174072.754290983080864)
Pixel Size = (4.777314267160000,-4.777314267160000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-6398696.511,-1174072.754) ( 57d28'49.69"W, 10d29'16.12"S)
Lower Left (-6398696.511,-1193640.634) ( 57d28'49.69"W, 10d39'38.19"S)
Upper Right (-6379128.632,-1174072.754) ( 57d18'16.87"W, 10d29'16.12"S)
Lower Right (-6379128.632,-1193640.634) ( 57d18'16.87"W, 10d39'38.19"S)
Center (-6388912.571,-1183856.694) ( 57d23'33.28"W, 10d34'27.20"S)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Blue
Overviews: 2048x2048, 1024x1024, 512x512
Band 2 Block=512x512 Type=UInt16, ColorInterp=Green
Overviews: 2048x2048, 1024x1024, 512x512
Band 3 Block=512x512 Type=UInt16, ColorInterp=Red
Overviews: 2048x2048, 1024x1024, 512x512
Band 4 Block=512x512 Type=UInt16, ColorInterp=Undefined
Overviews: 2048x2048, 1024x1024, 512x512
Band 5 Block=512x512 Type=UInt16, ColorInterp=Alpha
Overviews: 2048x2048, 1024x1024, 512x512
This command gives you some basic information into the metadata of your quads. For example, it tells you how many bands your quad has (4 bands + alpha) and in what projection it is in (WGS84). Now that we have inspected the tiff files with the above command, we will run the following command to stitch them together.
!gdal_merge.py -v data/L15-0697E-0963N.tif data/L15-0698E-0963N.tif -o output/merged.tif
The -v flag in the command allows us to see the output of the mosaicing operations as they are done. We can see in the verbose output above that the mosaicing operation is fairly simple: the utility script is basically copying a range of pixels for each band in each scene over to the designated output file.
Processing file 1 of 2, 0.000% completed in 0 minutes.
Filename: L15-0697E-0963N.tif
File Size: 4096x4096x5
Pixel Size: 4.777314 x -4.777314
UL:(-6398696.510914,-1174072.754291) LR:(-6379128.631675,-1193640.633529)
Copy 0,0,4096,4096 to 0,0,4096,4096.
Copy 0,0,4096,4096 to 0,0,4096,4096.
Copy 0,0,4096,4096 to 0,0,4096,4096.
Copy 0,0,4096,4096 to 0,0,4096,4096.
Copy 0,0,4096,4096 to 0,0,4096,4096.
Processing file 2 of 2, 50.000% completed in 0 minutes.
Filename: L15-0698E-0963N.tif
File Size: 4096x4096x5
Pixel Size: 4.777314 x -4.777314
UL:(-6379128.631675,-1174072.754291) LR:(-6359560.752437,-1193640.633529)
Copy 0,0,4096,4096 to 4096,0,4096,4096.
Copy 0,0,4096,4096 to 4096,0,4096,4096.
Copy 0,0,4096,4096 to 4096,0,4096,4096.
Copy 0,0,4096,4096 to 4096,0,4096,4096.
Copy 0,0,4096,4096 to 4096,0,4096,4096.
Let's run gdalinfo again to see the spatial metadata of the mosaiced file.
Driver: GTiff/GeoTIFF
Files: merged.tif
Size is 8192, 4096
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
Origin = (-6398696.510913709178567,-1174072.754290983080864)
Pixel Size = (4.777314267160000,-4.777314267160000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-6398696.511,-1174072.754) ( 57d28'49.69"W, 10d29'16.12"S)
Lower Left (-6398696.511,-1193640.634) ( 57d28'49.69"W, 10d39'38.19"S)
Upper Right (-6359560.752,-1174072.754) ( 57d 7'44.06"W, 10d29'16.12"S)
Lower Right (-6359560.752,-1193640.634) ( 57d 7'44.06"W, 10d39'38.19"S)
Center (-6379128.632,-1183856.694) ( 57d18'16.87"W, 10d34'27.20"S)
Band 1 Block=8192x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8192x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8192x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8192x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8192x1 Type=UInt16, ColorInterp=Undefined
As you can see above, the size of the merged raster has changed from 4096,4096 to 8192,4096.
2. Clip to an AOI¶
Now that we have succesfully mosaiced the quads, we will now clip the merged file to our exact area of interest. For this step, we will run the following command:
!gdalwarp -of GTiff -cutline data/clip_aoi.geojson -crop_to_cutline output/merged.tif output/merged_cropped.tif
3. Perform NDVI on the clipped mosaic¶
For the NDVI calculation, we will be taking a look at the 3rd and 4rth band of the merged file. Jut as a note: the band rendering for analytic quads are the following: Band 1: Blue Band 2: Green Band 3: Red Band 4: NIR Band 5: Alpha
from osgeo import gdal
import numpy as np
g = gdal.Open("merged.tif")
if g is None:
raise IOError("Couldn't open merged.tif")
b3 = g.GetRasterBand(3).ReadAsArray().astype(np.float32)
b4 = g.GetRasterBand(4).ReadAsArray().astype(np.float32)
ndvi = (b4 - b3) / (b4 + b3)
ndvi = ndvi.astype(np.float32)
drv = gdal.GetDriverByName("GTiff")
dst_ds = drv.Create("output_file.tif", g.RasterXSize, g.RasterYSize, 1, gdal.GDT_Float32, options=["COMPRESS=LZW"])
dst_ds.GetRasterBand(1).WriteArray(ndvi)
dst_ds = None
Good Job! You have finished the workflow!
Rate this guide: