Python Landsat Tutorial

Requirements

This tutorial is designed to show how to work with Landsat 8 files in WASDI. Details of the Landsat 8 mission and/or guidelines on how to configure your own environment are out of the scope of this tutorial In this tutorial we use PyCharm as a free Python Development tool, but the code can be executed on every different Python environment.

Note

This tutorial requires gdal working in your python env: we know this can be tricky.

For Windows 10, we suggest following this tutorial:

https://opensourceoptions.com/blog/how-to-install-gdal-for-python-with-pip-on-windows/

The key is that here:

https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

you can find different pre-build GDAL “images” that can be installed directly asking pip to use the file.

Overview

In this tutorial, we will learn how to work with Landsat 8 Images in WASDI. The tutorial will implement a processor that takes as input:

  • Name of a Landsat 8 file

  • List of Bands to extract, by default B5 and B4

  • Resolution to use

The processor will open the Landsat image, create a tiff with only the selected bands, reprojected to the desired resolution, and then it will calculate NDVI, assuming that the first extracted band is NIR and the second is RED.

The equation to compute NDVI is: [NIR-RED]/[NIR+RED].

Note

It is mandatory that at least 1 Landsat image is imported in the workspace (using the WASDI web interface, i.e. Search) BEFORE running this processor. Check Wasdi Web Platform access and basic usage for more general info on this.

Setup

Open PyCharm and start a new project.

../_images/openPyCharmAndStartANewProject.png

Name it “Landsat8Tutorial” (or however you wish, just remember to be coherent). You may wish to create a new virtual environment or use an existing one. Uncheck the option for creating a “main.py” welcome script (or, at least, remember to delete it later on).

Let’s install the library we need. In the terminal write:

pip install wasdi

and hit enter

../_images/pipInstallWasdi.png

Note

Remember you need also to have gdal installed. If you previously installed wasdi, you may wish to update it by adding the –upgrade flag, i.e.:

pip install --upgrade wasdi

Create first files

Now we need to create these three fundamental files (right click on the Projec Icon, new -> …):

  • myProcessor.py: create a python file, then call it myProcessor.py

  • config.json: create a file, then call it config.json (PyCharm will recognize automatically it is a JSON file)

  • params.json: create a file, then call it params.json (PyCharm will recognize automatically it is a JSON file)

../_images/createPythonFile.png

Create python file

../_images/createMyProcessor_py.png

Call it myProcessor.py

../_images/createFile.png

Create a json file

../_images/createConfig_json.png

Call it config.json

Create a json file

../_images/helloWASDIWorldDebug1.png

Call it params.json

Next, point your browser to wasdi.net, log in, go in the Workspaces Section and create a new workspace. Call it “Landsat8Tutorial”.

Go to the search section and select L8 data type and a bounding box in Europe

../_images/searchL8Image.jpg

Select one image and click on the + button to add the image to the Landsat8Tutorial Wokspace

../_images/addToWS.jpg

Come back to the edit section, and check that WASDI has been able to fetch the image.

../_images/imageInWorkspace.jpg

Take note of the file you imported, we will need it later. For this tutorial we assume:

LC08_L1GT_196029_20211227_20211227_01_RT

but this can be changed with any image you imported.

Leave the browser open on that page, we will need it later on.

First lines

Let’s begin by editing the config.json file. It is a JSON file, containing the user credentials and some fundamental parameters to get you started (see Wasdi Libraries Concepts):

{
  "USER": "your user name here",
  "PASSWORD": "your password here",
  "PARAMETERSFILEPATH": "./params.json"
  "WORKSPACE": "AdvancedTutorialTest"
}

Note

please, keep this file for yourself. You should never give this file to anyone else, and you do not need to upload to WASDI, as we’ll see later on. You just need this file in your project for working with the WASDI python library. Use this file to change the workspace where you want to work.

Let’s then edit params.json file. It is a JSON file that represents the inputs needed by our processor. The WASDI Developer can decide what parameters are needed; each parameter has a unique name within the processor. Each parameter can be of different types (i.e. Strings, Integers, Float, Arrays, Complex Objects…). params.json is where you declare and valorize your inputs. The same inputs will be avaiable in the WASDI Web Interface when publishing the processor.

{
  "BANDS": ["B5", "B4"],
   "RESOLUTION": "30",
   "L8FILE": "LC08_L1GT_196029_20211227_20211227_01_RT.zip"
}

Now, open myProcessor.py, create a main and a method called run. The latter is required for WASDI to work (more on that later on).

Note

These are two requirements necessary to use WASDI:

  • have a python file called myProcessor.py

  • have a function called run() (no params) within myProcessor.py

After that, you can include as many python files as you need, regardless their organization in directories. You just need to have a myProcessor.py with a method run() as entry point.

The main method will initiate the WASDI library and call the run method:

import wasdi


def run():
    pass


if __name__ == '__main__':
    wasdi.init("./config.json")
    run()

As you can see, we call wasdi.init and pass the relative path of the config file to it.

../_images/wasdi_init.png

Let’s debug to see the effects of this.

Note

If a file main.py was created automatically for you, remember to define another debug configuration. The easiest way to do so is by right clicking on your code and select Debug ‘myProcessor.py’.

../_images/helloWASDIWorldDebug0.png

If the setup is correct so far, we should see the output from the wasdi library that shows the initialization has gone well. Something like this:

[INFO] _loadParams: wasdi could not load param file. That is fine, you can still load it later, don't worry
[INFO] waspy.init: returned session is: 0d3f3ef1-f4c3-4202-9015-6ca17fc21cc7
[INFO] waspy.init: WASPY successfully initiated :-)
[INFO] waspy.printStatus: user: username@email.address
[INFO] waspy.printStatus: password: ***********
[INFO] waspy.printStatus: session id: 0d3f3ef1-f4c3-4202-9015-6ca17fc21cc7
[INFO] waspy.printStatus: active workspace: 4f541d2c-4b29-445b-9869-9c8d185932ce
[INFO] waspy.printStatus: workspace owner: username@email.address
[INFO] waspy.printStatus: parameters file path: [...]/params.json
[INFO] waspy.printStatus: base path: C:\Users\username\.wasdi\
[INFO] waspy.printStatus: download active: True
[INFO] waspy.printStatus: upload active: True
[INFO] waspy.printStatus: verbose: True
[INFO] waspy.printStatus: param dict: {'BANDS': ['B5', 'B4'], 'RESOLUTION': '30', 'L8FILE': 'LC08_L1GT_196029_20211227_20211227_01_RT.zip'}
[INFO] waspy.printStatus: proc id:
[INFO] waspy.printStatus: base url: http://www.wasdi.net/wasdiwebserver/rest
[INFO] waspy.printStatus: is on server: False
[INFO] waspy.printStatus: workspace base url: http://www.wasdi.net/wasdiwebserver/rest
[INFO] waspy.printStatus: session is valid :-)

If you have the same situation, we are configured and ready to start!!

Extract Bands

The first step of our processor will be to extract the bands from the L8 image. WASDI ingest L8 images as a .zip file. Each .zip file contains different .tif images, one for each band, and some other files. We want to implement a function that takes as input the name of the L8 zip file, a list of bands, a resolution and that then creates a new .tif file with only the extracted bands at the desired resolution. The L8 bands are:

  • B1 - Coastal aerosol 30m

  • B2 - Blue 30m

  • B3 - Green 30m

  • B4 - Red 30m

  • B5 - Near Infrared (NIR) 30m

  • B6 - SWIR 1 30m

  • B7 - SWIR 2 30m

  • B8 - Panchromatic 15m

  • B9 - Cirrus 30m

  • B10 - Thermal Infrared (TIRS) 1 100m

  • B11 - Thermal Infrared (TIRS) 2 100m

Our function is implemented like this:

def extractBands(sFile, asBands, sResolution="30"):
   """
   Extracts some bands from the L8 zip file into a multiband tiff file at the specified resolution
    Bands are
    B1 - Coastal aerosol 30m
    B2 - Blue  30m
    B3 - Green 30m
    B4 - Red   30m
    B5 - Near Infrared (NIR) 30m
    B6 - SWIR 1 30m
    B7 - SWIR 2 30m
    B8 - Panchromatic 15m
    B9 - Cirrus 30m
    B10 - Thermal Infrared (TIRS) 1 100m
    B11 - Thermal Infrared (TIRS) 2 100m

    :param sFile: name of the Landsat 8 file
    :param asBands: array of string with the names of the bands to extract
    :param sResolution: resolution as a string is in meteres
    :return Returns the name of the new tiff file
   """

   # Output File Name that will be returned
   sOutputTiffFile = ""

   try:
       # Prepare the name a .vrt file that will be used to extract bands from the zip
       sOutputVrtFile = sFile.replace(".zip", ".vrt")
       # Prepare the name of the ouptut tif file
       sOutputTiffFile = sFile.replace(".zip", ".tif")

       # Get the Local Path of the input Landsat file
       sLocalFilePath = wasdi.getPath(sFile)

       # Get the path of the output files
       sOutputVrtPath = wasdi.getPath(sOutputVrtFile)
       sOutputTiffPath = wasdi.getPath(sOutputTiffFile)

       # Prepare an array of bands called BXX.TIF
       asBandsTiff = [b + '.TIF' for b in asBands]

       # Open the zip file
       with zipfile.ZipFile(sLocalFilePath, 'r') as zf:
           # Get all the files in the zip
           asZipNameList = zf.namelist()
           # Take from the files in the zip, the ones that match the BXX.TIF naming schema we are searching
           asBandsL8 = [name for name in asZipNameList for band in asBandsTiff if band in name]

           # Create the zip path of the files we want to extract
           asBandsZip = ['/vsizip/' + sLocalFilePath + '/' + band for band in asBandsL8]

           # Create an array that has the names of the files to extract in the order required by the asBands array in input
           asOrderedZipBands = []

           for sBand in asBands:
               for sZipBand in asBandsZip:
                   if sBand in sZipBand:
                       asOrderedZipBands.append(sZipBand)
                       break

           # Let gdal build a virtual file with our bands
           gdal.BuildVRT(sOutputVrtPath, asOrderedZipBands, separate=True)

           # Convert the vrt in tif with option  -tr sResolution sResolution to have all bands at the same res (ie -tr 30 30 to have at 30 meters)
           gdal.Translate(sOutputTiffPath, sOutputVrtPath, options="-tr " + sResolution + " " + sResolution)

           # we can remove the vrt file
           os.remove(sOutputVrtPath)
   except Exception as oEx:
       wasdi.wasdiLog("extractBands EXCEPTION")
       wasdi.wasdiLog(repr(oEx))
       wasdi.wasdiLog(traceback.format_exc())
   except:
       wasdi.wasdiLog("extractBands generic EXCEPTION")

   # Return the output file name
   return sOutputTiffFile

Compute NDVI

The second step is to compute the NDVI starting for our extracted Tif file. To compute NDVI we need to access the NIR and RED bands and compute the formula: NDVI = NIR-RED/NIR+RED

def computeNDVI(sTiffFile, sNDVIOutputFile):
    """
    Compute ndvi assuming that in sTiffPath there is as band 1 NIR and band 2 RED
    :param sTiffFile: name of the input tiff file
    :param sNDVIOutputFile: name of the ouput file with ndvi
    :return: full path of sNDVIOutputFile
    """

    # Open the tiff file: we assume it has two bands
    oDataset = gdal.Open(wasdi.getPath(sTiffFile))

    if not oDataset:
        wasdi.wasdiLog("Impossible to get Dataset from " + sTiffFile)
        return ""

    # Get the dimension of the bands in input
    [iCols, iRows] = oDataset.GetRasterBand(1).ReadAsArray().shape
    # Create gdal GeoTiff driver
    oDriver = gdal.GetDriverByName("GTiff")
    # Create a new Ouput file, same dimension of the input, compressed and with type float32.
    oOutDataFile = oDriver.Create(wasdi.getPath(sNDVIOutputFile), iRows, iCols, 1, gdal.GDT_Float32, ['COMPRESS=LZW', 'BIGTIFF=YES'])

    # set to the output same geotransform as input
    oOutDataFile.SetGeoTransform(oDataset.GetGeoTransform())
    # set to the output same projection as input
    oOutDataFile.SetProjection(oDataset.GetProjection())

    # We assume NIR = band1, RED = band2
    oNIR = oDataset.GetRasterBand(1)
    oRED = oDataset.GetRasterBand(2)

    # Convert the band values in a numpy array
    adNIRBandArray = numpy.array(oNIR.ReadAsArray())
    adREDBandArray = numpy.array(oRED.ReadAsArray())
    # Force data to be float
    adNIRBandArray = adNIRBandArray.astype(float)
    adREDBandArray = adREDBandArray.astype(float)
    # Compute NDVI formula, where is not nan
    adNDVIBandArray = numpy.where((adNIRBandArray + adREDBandArray!=0), (adNIRBandArray-adREDBandArray)/(adNIRBandArray+adREDBandArray), 0)

    # Write the new calulated NDVI to ouput file band 1
    oOutDataFile.GetRasterBand(1).WriteArray(adNDVIBandArray)
    # We assume 0 as no data
    oOutDataFile.GetRasterBand(1).SetNoDataValue(0)

    # saves to disk!!
    oOutDataFile.FlushCache()
    wasdi.wasdiLog("Saved " + sNDVIOutputFile)

    # Clean memory
    oNIR = None
    oRED = None

    # Return the name of our NDVI create file
    return sNDVIOutputFile

This tutorial shows an NDVI as a sample, but is clear that with this technique you can manipulate L8 data to fit your needs.

Main Function

Now the main operations are ready, we just need to put it all togheter.

def run():
    wasdi.wasdiLog("Landsat tutorial v.1.0")

    # Read from params the bands we want to extract and the resolution
    asBands = wasdi.getParameter("BANDS", ["B5", "B4"])
    sResolution = wasdi.getParameter("RESOLUTION", "30")
    sL8File = wasdi.getParameter("L8FILE", "LC08_L1GT_196029_20211227_20211227_01_RT.zip")

    # Call extract bands
    sTiffFile = extractBands(sL8File, asBands, sResolution)

    # Prepare the output NDVI name
    sNDVIFile = sTiffFile.replace(".tif", "_NDVI.tif")

    # Call compute NDVI
    computeNDVI(sTiffFile, sNDVIFile)

    # Add the file to the WASDI workspace
    wasdi.addFileToWASDI(sNDVIFile, "NDVI")

You can now test your processor. Remember that, at the first time you will debug it locally, WASDI will take some time to download for you the L8 file you are using. All is done automatically and only once, when needed.

In the same way, when you add the file to WASDI, the lib will updload for your result to the cloud:

[INFO] waspy._internalAddFileToWASDI( LC08_L1GT_196029_20211227_20211227_01_RT_NDVI.tif, False )
[INFO] waspy._internalAddFileToWASDI: remote file is missing, uploading
upload LC08_L1GT_196029_20211227_20211227_01_RT_NDVI.tif
uploadFile: uploading file to wasdi...
uploadFile: upload complete :-)
[INFO] waspy._internalAddFileToWASDI: file uploaded, keep on working!
[INFO] Running Locally, will not update status on server

Now that the core of our processor is done, lets make it a little bit more WASDI-integrated. We want to give some feedback to the user while the app is runnig and we do this using:

  • wasdi.wasdiLog: locally just a print to console, when on the server, it sends the logs to the web user interface

  • wasdi.updateProgressPerc: when on the server, updates the progress bar of the processor

  • wasdi.setPayload: allows to save a user-defined object associated to the processor run

def run():
    wasdi.wasdiLog("Landsat tutorial v.1.0")

    # Read from params the bands we want to extract and the resolution
    asBands = wasdi.getParameter("BANDS", ["B5", "B4"])
    sResolution = wasdi.getParameter("RESOLUTION", "30")
    sL8File = wasdi.getParameter("L8FILE", "LC08_L1GT_196029_20211227_20211227_01_RT.zip")

    wasdi.wasdiLog("Calling extract bands")
    # Call extract bands
    sTiffFile = extractBands(sL8File, asBands, sResolution)

    wasdi.updateProgressPerc(30)
    wasdi.wasdiLog("Calculating NDVI")

    # Prepare the output NDVI name
    sNDVIFile = sTiffFile.replace(".tif", "_NDVI.tif")

    # Call compute NDVI
    computeNDVI(sTiffFile, sNDVIFile)
    wasdi.updateProgressPerc(80)

    wasdi.wasdiLog("Adding " + sNDVIFile + " to the workspace")
    # Add the file to the WASDI workspace
    wasdi.addFileToWASDI(sNDVIFile, "NDVI")

    # Create the payload object
    aoPayload = {}
    # Save the inputs that we received
    aoPayload["inputs"] = wasdi.getParametersDict()
    # Save the output we created
    aoPayload["output"] = sNDVIFile
    # Save the payload
    wasdi.setPayload(aoPayload)

    # Close the process setting the status to DONE
    wasdi.updateStatus("DONE", 100)

Welcome to Space, Have fun!