July 26 2021

Beginners QGIS Workflow with NICFI data

Introduction

This tutorial will cover basic raster operations in QGIS. You will first mosaic two basemap quads and then clip the mosaic to an area of interest. After clip, you are going to use Raster Calcualtor to perform NIDVI. The steps in this workflow are

  1. Mosaic basemap quads
  2. Clip to an AOI
  3. Create NDVI Raster

You can choose to follow all three of the steps, or only follow one step. The workflow is meant for beginners starting to work with NICFI data.

Before you start

  • You will need QGIS installed on your machine to perform the raster operations in this tutorial. This tutorial is using QGIS version 3.12 but the raster operations should be available in all QGIS versions.
  • You will need to download at least two basemap quads from Planet.
  • You will need an area of interest that you want to clip the mosaic. Area of interest can be in GeoJSON or shapefile, and should work with all QGIS supported vector formats.

While this is a simple tutorial, it is recommended that you have basic familiarity with GIS concepts. You can find more on working with raster data in QGIS here.

1. Mosaic Basemap Quads

To start, open a QGIS project and add your GeoTiff basemap quads (It is assumed that you have the two basemap quads stored locally in your machine).

Loading Raster Data:

Layer > Add Layer > Add Raster Layer

On the pop up window, under source, click the "..." and choose your GeoTiff quads. More than one layer can be loaded at the same time by holding down the Ctrl and Shift key and clicking on multiple items. Click Add to add the quads to your project.

1.1 Merge two quads together:

  1. Click on Raster
  2. Hover over Miscellaneous
  3. Click on Merge
  4. On the Merge dialogue
    1. In the Input Layers click on ...
    2. In the Output data type choose Int16
    3. Click Run

merge.gif

Good Job! You have merged multpile quads together!

1.2 Save the Merge output:

  1. Right click on the Merge output
  2. Export
  3. Click Save As
  4. Navigate to where you want to save your Merged GeoTiff and provide name
  5. Choose CRS (Planet mosaic quads are ESPG:3857 - WGS 84 / Pseudo)
  6. Click OK

2. Clip to an AOI

In this step, you are are going to clip the merged GeoTIFF to a area of interst (e.g. shapefile or GeoJSON).

Loading Shapefile Data:

Layer > Add Layer > Add Vector Layer

On the pop up window, under source, click the "..." and choose your clip file. As with the GeoTIFFs, more than one vector layer can be loaded at the same time by holding down the Ctrl and Shift key and clicking on multiple items. Click Add to add the quads to your project.

Example of merged quads before clip.

before_clip.png

2.1 Clip the merged raster to features

  1. Right click on Raster
  2. Hover over Extraction
  3. Click on Clip Raster by Mask Layer
  4. On the Clip Raster By Extent dialogue
    1. Select Input Layer (merged_quads.tif)
    2. Selet Mask Layer (clip_features.shp)
    3. Click Run

clip.gif

2.2 Save the Merge output:

  1. Right click on the Clipped (mask) output
  2. Export
  3. Click Save As
  4. Navigate to where you want to save your Merged GeoTiff and provide name
  5. Choose CRS (Planet mosaic quads are ESPG:3857 - WGS 84 / Pseudo)
  6. Click OK

Example of merged quads after clip.

after_clip.png

3. Create NDVI Raster

Finally, you will calculate NDVI (Normalized Difference Vegetation Index) and create a new raster with these values.

3.1 NDVI Overview

NDVI, developed by a NASA scientist named Compton Tucker in 1977, is commonly used to assess whether an area contains live green vegetation or not. It can show the difference between water and plants, bare soil and grass, whether plants are under stress, and what lifecycle stage a crop is in.

It compares how much more near-infrared light is reflected by chlorophyll vs visible red light and is given as the following equation:

ndvi_equation.png

3.2 Calculating NDVI in QGIS

  1. Click on Raster
  2. Choose Raster Calculator
  3. In the Raster Calculator Expression
    1. Type ( "clipped_quads@4" - "clipped_quads@3" ) / ( "clipped_quads@4" + "clipped_quads@3" ) where clipped_quads is the raster image you want to perfrom NDVI.
  4. Save your NDVI by choosing the path in Output layer

ndvi_calc.png

Nice work! The final step is to make the color ramp vegetation friendly.

3.3 Edit Color Ramp

  1. Right click on the NDVI layer
  2. Choose Properties
  3. On the left side click Symbology
  4. In the Render Type choose singleband pseudocolor
  5. For mode pick equal interval
  6. Click Apply and then OK

color_ndvi.gif

Good work! You have finished the whole workflow!

July 26 2021

Introduction

This tutorial will cover basic raster operations in ArcGIS Pro. In this workflow you will:

  1. Mosaic two basemap quads
  2. Clip the mosaic to a polygon feature class
  3. Create an NDVI raster

Before you start

You will need ESRI ArcGIS Pro or ArcMap to perform the raster operations demonstrated in this tutorial. Additionally, you'll need at least an ArcGIS Pro Basic lisence to use the necessary tools. More information on ESRI license requirements can be found here.

You will also need to have previously downloaded at least two basemap quads from Planet. Information on downloading basemap quads from Planet's Basemaps Viewer web application can be found here. Finally, you will also need a polygon feature class that will be used to clip the mosaic you create.

While this is a simple tutorial, it is recommended that you have basic familiarity with GIS concepts.

1. Mosaic Basemap Quads

After downloading the two basemaps quads and creating a new project in ArcGIS Pro, click the Add Data button on the Map ribbon to add the basemap quads and clip feature class to your project.

add_data.png

load_data.png

Next, navigate to the Analysis ribbon and go to Raster Functions. ESRI has extensive documentation on the Raster Functions toolset if you are interested in learning more.

raster_fn.png

Under Data Management click Mosaic Rasters to open the tool dialog box. Alternatively, you can find the tool using the search bar atop the raster Raster Functions pane.

data_mngmt.png

With the Mosaic Rasters tool dialog open, select your basemap quads for the Raster parameter. Under Operation, select Max. To learn more about any tool in ArcGIS Pro, click the question mark in the upper right-hand corner of the tool dialog.

mosaic_rasters.png

After clicking Create new layer at the bottom of the tool dialog box, your new mosaic image will appear in the map view.

2. Clip Mosaic

Next, you will clip your newly created mosaic to the polygon feature class.

From the Raster Functions pane select the Clip tool that is also located under Data Management.

In the following tool dialog box, select the mosaic for the Raster parameter, a Clipping Type of "Outside" (this will remove pixels outside of the polygon clip feature class), and select the polygon feature class for the Clipping Geometry parameter.

Be sure to check the box next to Use input features for clipping geometry, which will insure that the mosaic is clipped to the input feature class.

clip_dialog.png

After clicking Create new layer at the bottom of the tool dialog box, your newly clipped mosaic will appear in the Map view.

clip_mosaic.png

Nice work! You have successfully mosaiced two basemap quads and clipped the mosaic.

3. Create NDVI Raster

Finally, you will calculate NDVI (Normalized Difference Vegetation Index) and create a new raster with these values.

I. NDVI Overview

NDVI, developed by a NASA scientist named Compton Tucker in 1977, is commonly used to assess whether an area contains live green vegetation or not. It can show the difference between water and plants, bare soil and grass, whether plants are under stress, and what lifecycle stage a crop is in.

It compares how much more near-infrared light is reflected by chlorophyll vs visible red light and is given as the following equation:

ndvi_equation.png

II. Calculating NDVI in ArcGIS Pro

There are a number of ways to calculate NDVI in ArcGIS Pro. In this tutorial, you'll be utilizing the native NDVI tool which is located in the Raster Functions pane.

In the following tool dialog box, select the clipped image for the Raster parameter.

Under Visible Band ID, select 3. Under Infrared Band ID, select 4 For this basemap, band 3 is the red band, and band 4 is the near-infrared band.

Optionally, check the box next to Scientific Output, which will produce NDVI values between -1 and 1.

NDVI_dialog.png

After clicking Create new layer at the bottom of the tool dialog box, your new NDVI raster will appear in the Map view. You may wish to adjust the symbology of the raster.

final.png

Nice work! You have successfully created an NDVI raster from the clipped mosaic!

July 26 2021

Introduction

In this tutorial, you will do the following:

  1. Mosaic Basemap Quads
  2. Clip to an AOI
  3. 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!