Satellite imagery brought us the capacity to see the land surface on recent years but we haven’t been so successful to understand land cover dynamics and the interaction with economical, sociological and political factors. Some deficiencies were found on the use of GIS commercial software, but there are other limitations in the way we apply logical and mathematical processes to a set of satellite imagery. Working with geospatial data on Python gives us the capability to filter, calculate, clip, loop, and export raster or vector datasets with an efficient use of the computational power providing a bigger scope on data analysis.

This tutorial shows the complete procedure to create a land cover change raster from a comparison of generated vegetation index (NDVI) rasters by the use of Python and the Numpy and GDAL libraries. Contours of land cover change where generated with some tools of GDAL and Osgeo and an analysis of deforestation were done based on the output data and historical images from Google Earth.

Current GIS desktop applications are fully capable of this spatial management and analysis when the amount of raster images is limited; however when we deal with high amount of images the spatial processing on a graphical user interfase (GUI) can be slow and most commonly impractical. The use of programming / processing languages like Python and advanced spatial libraries as GDAL (gdal.org) helps on the spatial data transformation on a more abstract and effective way. This tutorial shows the complete procedure to clip the complete set of bands from a Landsat 8 image and store them with a suffix on every band file on another folder.

The tutorial is done on a interactive Python programming platform called Jupyter Notebook. The input files: raster bands and area of interest (AOI) shapefile need to be on the same system of reference (SRC), otherwise the GDAL library cannot locate the spatial data on the right position. The tutorial shows the procedure for the whole set of band form a Landsat 8 image, an example for a single band is provided on the scripts of the input data. Finally the tutorial shows the complete and clipped raster on a GIS desktop software as QGIS.

Satellite images are georasters, these images are a regular array of columns and rows (a matrix per band) with a georeferenciation. Python is a programming and data analysis language very versatile for the matrix algebra with the Numpy library, however there was no efective and simple way to process a georaster until the development of the Rasterio package.

Rasterio is a library to open, write, explore and analyze georasters in Python. The library uses GeoTIFF images along with other formats and is capable to work with satellite images, digital elevation models, and drone generated imagery.

This tutorial show the complete procedure to analyse the NDVI from a Landsat 8 image with Python 3 and Rasterio. The scripting and representation was performed on a interactive enviroment called Jupyter Notebook, finally the result georaster was opened in QGIS and compared with some background images.

Rasterio is a Python library that allows to read, inspect, visualize and write geospatial raster data. The library uses GeoTIFF and other spatial raster formats and is capable of working with satellite imagery, digital elevation models, and drone imagery data products. Rasterio allows you to import a single band or multiband geospatial raster in a interactive Python enviroment as Jupyter notebook, the library can keep the “duality” of the geospatial raster, that means, it can handle the location and resolution parameters as well as the matrix values of the gridded elements.

This tutorial shows some basic procedures to explore a multiband Sentinel 2 granule with Python 3 and Rasterio on a Jupyter Notebook. The tutorial shows the commands to identify the raster array dimensions and the geospatial referencing parameters, make representation of each visible band and export band composites as true color and false color geoespatial rasters in Tiff format.

Digital elevation models (DEMs) from satellite interpretation (Aster DEM or Alos Palsar) come with “sinks” from errors in the elevation interpretation, raster resolution or reprojection. There is a need to correct those rasters in order to interpret the hydrological features. This tutorial show the process to condition a digital elevation model (DEM) dowloaded from a NASA/USGS server (gdex.cr.usgs.gov) with the Pysheds library of Python. The tutorial was done on a Jupyter Notebook, input files and scripts are attached on the final part of the post.

Calculate volumes is an easy task on QGIS 3 with the use of SAGA GIS tools under the Processing toolbox. Raster volumes can be calculated above or below a base level and two other specific calculation, input rasters need to be in UTM. These type of calculation is useful for earth movement, reservoir design, risk analysis and other spatial analysis.

The tutorial presented on this post shows the complete procedure to calculate the raster volume above a base level on QGIS3.

What would happen if we shift our GIS geoprocessing to Python? What would happen if we treat our raster and vector spatial data as objects and variables on a Python 3 script? Then we can ask ourselves if it is neccessary to reinvent the wheel, it is necessary to change a workflow that already work on a GIS software.

There is a simple answer to this dilemma: More control

Working with Python give us more control on the geoprocessing itself since we leave the Graphical User Interface (GUI) with its icons, buttons and dialog boxes. With Python running on a Jupyter Notebook, we can link with specific files, define geoprocess and it options, make plots of draft and final data, and export results to vector/raster SIG formats. There are other advantages of spatial analysis in Python which are the reproducibility and the processing speed.

In GIS the objects are related to a spatial location. Usually there is no need to modify a location since spatial data comes from field work or other surveys. When we work with CAD files as DXF (Autocad Drawing Exchange Format) files, sometimes the spatial data is available as a layout view and not as a model view. The layout data is on local coordinates, and at specific scale therefore we need to scale and translate the spatial data to “return” it to its original spatial location and extension.

We tried to make this job with Inkscape and QGIS3, but we were unsuccessful to complete the scale and translation. While working with the DXF in Inkscape, each object selection, layer order combination, object ungrouping took several minutes and the results were poor. At the lowest motivational stage of this spatial request, we thought that it might me something in Python that can be useful for this.

On a internet search, we found that the spatial version of the Pandas data analysis library: Geopandas was capable not only to open the DXF files, but also to scale, translate, and filter spatial data according to specific criteria. Geopandas is capable to export spatial data in different formats and to plot data interactively on a Jupyter Notebook.

This tutorial shows the procedure to open a DXF file in Python pandas, perform scale and translation to place the spatial features on their original position, filter unwated objects on the layout view and export results to QGIS3 as shapefile.