Articles

QGIS Model Designer offers powerful functionality for automating workflows. Many workflows are based on conditional logic that can generate different outputs based on the inputs provided. This post shows how to setup a QGIS model that can skip certain steps if those optional inputs are not provided. This is achieved using the Conditional branch algorithm in the QGIS Model Designer. Here’s the example workflow we want to implement Given a buildings layer and a boundary layer, create a layer with building polygons at least 50m away from the boundary. This will be the final output if no other inputs are provided. Optionally, if a roads layer has been provided, further process the data and create a layer with buildings polygons that are both 50m away from the boundary and within 25m of a road segment. The key is to use an expression that checks if the roads input is NULL and selects the appropriate branch of the model to process. Below is the diagram showing the final model. We will go through the step-by-step process of building this model. You will find a link to download a GeoPackage file containing the final model and the sample dataset at the end of the post. …

Welcome to #30DaysOfQGIS! We have launched our Advanced QGIS course on YouTube and have designed this challenge to help you master QGIS! Spend 30 minutes each day for the next 30 days to level-up your QGIS skills. This course is the result of my 15+ years of experience using QGIS for large-scale spatial analysis and automating workflows. I am really excited to share this content with you - completely free. The material is full of tips, tricks and challenges that will make your learning journey fun and rewarding! All you have to do is show up everyday and spend half an hour watching the videos and practicing the exercise. Ready for #30DaysOfQGIS? Read on to know the details. …

When automating GIS workflows, one often needs to automate the creation of cartographic outputs. The QGIS Model Designer allows you to build a workflow by combining multiple Processing algorithms. QGIS now includes several algorithms under the Cartography category that allow you to integrate the map creation process within your model. In this post, we will explore the Export print layout as image (or Export print layout as PDF) algorithm to automate the creation of a fire map. We will build a model that will automatically Download the latest shapefile of active fires from FIRMS. Extract fires intersecting the continental US. Style the layer using a pre-configured QGIS style file. Render a pre-configured Print Layout. Whenever the model is run, it will output a map such as shown below. …

An important concept in spatial statistics is pixel weights. When calculating pixel statistics with a polygon, partial pixel overlaps are treated differently by different packages and you need to understand this to evaluate the accuracy of your results. Consider the following image. What is the correct answer? …

Modern versions of QGIS comes with a handy command-line utility called qgis_process. This allows you to access and run any Processing Tool, Script or Model from the Processing Toolbox on a terminal. This is very useful for automation since it doesn’t require you to open QGIS or manually click buttons. You can run the algorithms in a headless-mode and even schedule them to run them at specific times. This post covers the following topics How to launch qgis_process command on Windows, Mac and Linux. How to find the parameters and values for each algorithm and build your command Example showing how to do a spatial join on the command-line using the Join Attributes by Location algorithm Example showing how to run a model on the command line to automate a complex workflow Want to follow along? You can download the data package containing all the datasets used in this post. Before running each command, make sure to replace the paths in the commands with the paths on your computer. …

In this post, I describe how we can use built-in QGIS processing tools to create a workflow to split polygons into equal parts. Using a clever algorithm and Feature Iterator tool in the Processing Framework, we can easily split all features in a given polygon layer into equal parts. This workflow creates irregular polygons. There are many plugins that allow more control and can divide polygons into regularly sized equal sized polygons. See the complete list at the end of the post. The algorithm for splitting any polygon shape into equal parts is described in this post PostGIS Polygon Splitting by Paul Ramsey. We will see how this can be implemented in QGIS. …

In this post, I will outline techniques for computing weighted-centroids in both QGIS and Google Earth Engine. For a polygon feature, the centroid is the geometric center. It can also be thought of as the average coordinate of all points within the polygon. There are some uses cases where you may want to compute a weighted-centroid where some parts of the polygon gets higher ‘weight’ than others. The main use-case is to calculate a population-weighted centroid. One can also use Night Lights data as a proxy for urbanized population and calculate a nightlights-weighted centroid. Some applications include: Regional Planning: Locate the population-weighted centroid to know the most accessible location from the region. Network Analysis: For generating demand points in location-allocation analysis, you need to convert demands from regions to points. It preferable to compute populated-weighted centroids for a more accurate analysis. Do check out this twitter-thread by Raj Bhagat P for more discussion on weighted centroids. Different Weighted Centroids for the State of Karnataka, India (2015) …

I recently taught a 1-month long course on GIS Applications in Urban and Regional Planning. We explored how GIS can be applied to solve problems in 6 different thematic areas. In this post, I will outline different applications and show concrete examples of using open-datasets and open-source GIS software QGIS. Update! The full course material - including data packages and PDF handouts - is now available for free download. Download [qgis_urban_planning.zip] containing detailed step-by-step instructions and datasets. You can also access updated versions of some tutorials at qgistutorials.com. Here are the 6 thematic areas Land Use Planning and Management Crime Mapping and Analysis Solid Waste Management Urban Infrastructure and Utilities Urban Transportation Spatial Planning …

Rainfall is arguably the most frequently measured hydro-meteorological variable. It is a required input for many hydrological applications like runoff computations, flood forecasting as well as engineering design of structures. However, rainfall data in its raw form contain many gaps and inconsistent values. Therefore it is important to do rigorous validation of rain-gauge observation before incorporating them into analysis. World Bank’s National Hydrology Project (NHP) prescribes a set of primary and secondary validation methods in theManual of Rainfall Data Validation. Of particular interest to me are the spatial methods aimed to identify suspect values by comparison with neighboring stations. This spatial homogeneity test requires complex spatial and statistical data processing that can be quite challenging. I got an opportunity to work on a project that required automating the entire process of identifying and testing suspect stations. I ended up implementing it in QGIS using just Expressions and Processing Modeler. The whole solution required no custom code and was easily usable by an analyst in the QGIS environment. In this post, I will explain the details of the test and show you how you can use similar techniques for your own analysis. This workflow was presented as a live session on QGIS Open Day. You can watch the recording to understand the concepts and implementation. …

When working with raster data, you may sometimes need to deal with data gaps. These could be the result of sensor malfunction, processing errors or data corruption. Below is an example of data gap (i.e. no data values) in aerial imagery. Source Image: © Commission for Lands (COLA) ; Revolutionary Government of Zanzibar (RGoZ), Downloaded from OpenAerialMap. (Note: The data gap is simulated using a python script and is not part of the original dataset) …

In a previous post, I showed how to use the aggregate function to find neighbor polygons using QGIS. Using aggregate functions on the same layer allows us to easily do geoprocessing operations between features of a layer. This is very useful in many analysis that would typically require writing custom python scripts. Here I demonstrate another powerful function array_foreach that allows one to iterate over other features in QGIS expressions - enabling even more powerful analysis by writing just a single expression. …

Generating pseudo-random data is important for many aspects of research work. QGIS provides for many methods of generating random points to facilitate this. Recently, I ran into a problem where I wanted to generate random points inside a polygon - but I wanted the random points to have a certain distribution. I wanted to generate a dataset showing employee home locations for a company. Given a city boundary and the location of office, I wanted to have a point layer that showed where the employees lived. A simple ‘Random points within Polygon’ algorithm would not work here, since the distribution of points would not be uniform within the city. …

Table Joins are a way to join 2 separate layers based on a common attribute value. QGIS has a Join Attributes By Field Value algorithm that allows you to perform table joins. A limitation of this algorithm is that the field values must match exactly. If the values differ slightly - the join will fail. There are many times where you are trying to join 2 layers from different sources and they contain values which are similar but may not match exactly. Fortunately QGIS now has built-in fuzzy string matching functions that can be used - along with Aggregatefunction - to do table join based on fuzzy matches. …

Read my previous posts Summary Aggregate and Spatial Filters and Advanced Aggregate Expressions to Automate QA to learn more about the powerful aggregate function. The aggregate function in QGIS was designed to work with 2 separate input vector layers, but we can also make it work with a single layer. Essentially performing spatial queries for features within the layer. The following video explains the workflow. …

This post is the continuation of Summary Aggregate and Spatial Filters in QGIS. I have been exploring aggregate functions more and have found interesting ways to automate tasks in QGIS. One such example is helping automatically keeping track of feature edits to help with Quality Assurance (QA). …

QGIS expression engine has a powerful a summary aggregate function that can do spatial joins on the fly. This enables some very interesting uses. …

You may have seen a map where source and destination points are connected via curved lines. It is possible to create such a map in QGIS with a simple trick - using custom projections and densification of lines. I will outline the steps to create such a map. …

Rainfall is arguably the most frequently measured hydro-meteorological variable. It is a required input for many hydrological applications like runoff computations, flood forecasting as well as engineering design of structures. However, rainfall data in its raw form contain many gaps and inconsistent values. Therefore it is important to do rigorous validation of rain-gauge observation before incorporating them into analysis. World Bank’s National Hydrology Project (NHP) prescribes a set of primary and secondary validation methods in theManual of Rainfall Data Validation. Of particular interest to me are the spatial methods aimed to identify suspect values by comparison with neighboring stations. This spatial homogeneity test requires complex spatial and statistical data processing that can be quite challenging. I got an opportunity to work on a project that required automating the entire process of identifying and testing suspect stations. I ended up implementing it in QGIS using just Expressions and Processing Modeler. The whole solution required no custom code and was easily usable by an analyst in the QGIS environment. In this post, I will explain the details of the test and show you how you can use similar techniques for your own analysis. This workflow was presented as a live session on QGIS Open Day. You can watch the recording to understand the concepts and implementation. …

In a previous post, I showed how to use the aggregate function to find neighbor polygons using QGIS. Using aggregate functions on the same layer allows us to easily do geoprocessing operations between features of a layer. This is very useful in many analysis that would typically require writing custom python scripts. Here I demonstrate another powerful function array_foreach that allows one to iterate over other features in QGIS expressions - enabling even more powerful analysis by writing just a single expression. …

Table Joins are a way to join 2 separate layers based on a common attribute value. QGIS has a Join Attributes By Field Value algorithm that allows you to perform table joins. A limitation of this algorithm is that the field values must match exactly. If the values differ slightly - the join will fail. There are many times where you are trying to join 2 layers from different sources and they contain values which are similar but may not match exactly. Fortunately QGIS now has built-in fuzzy string matching functions that can be used - along with Aggregatefunction - to do table join based on fuzzy matches. …

Read my previous posts Summary Aggregate and Spatial Filters and Advanced Aggregate Expressions to Automate QA to learn more about the powerful aggregate function. The aggregate function in QGIS was designed to work with 2 separate input vector layers, but we can also make it work with a single layer. Essentially performing spatial queries for features within the layer. The following video explains the workflow. …

This post is the continuation of Summary Aggregate and Spatial Filters in QGIS. I have been exploring aggregate functions more and have found interesting ways to automate tasks in QGIS. One such example is helping automatically keeping track of feature edits to help with Quality Assurance (QA). …

QGIS expression engine has a powerful a summary aggregate function that can do spatial joins on the fly. This enables some very interesting uses. …

QGIS comes bundled with a simplified version of the Natural Earth Countries shapefile that is suitable for quick map-making. The layer can be loaded into your canvas by typing the keyword world in the coordinates bar. While this is useful, there is no single political map of the world that is accepted by every country of the world. There are many disputed international boundaries, and each country has its own version of accepted international boundaries. To allow mapmakers to adhere to local mapping regulations, Natural Earth also publishes Countries point-of-views shapefiles for many countries that depict the world map according to each country’s law and/or local conventions. We provide a simple script to replace the bundled world map with your country’s point-of-view layer. …

Welcome to #PyQGISChallenge - Master QGIS Python Development in 30 Days! We are launching our PyQGIS Masterclass course on YouTube and have designed this challenge to help you learn how to customize QGIS using Python with scripts, custom algorithms, actions and plugins! Spend 30 minutes each day for the next 30 days to level-up your QGIS skills. This course is the result of my 15+ years of experience doing QGIS development - including building enterprise-grade plugins and deploying QGIS to thousands of users. I am really excited to share this content with you - completely free. We will be posting short videos everyday and cover the full course material step by step. The material is designed to help you slowly ramp up and learn complex concepts! All you have to do is show up everyday and spend half an hour watching the videos and practicing the exercises. At the end, you can take up a mini-project and apply your newly acquired skills. Ready for #PyQGISChallenge? Read on to know the details. This is an advanced course that assumes good working knowledge of both Python and QGIS. If you are new to programming, complete our Python Foundation for Spatial Analysis course first. …

K-Means Clustering is a popular algorithm for automatically grouping points into natural clusters. QGIS comes with a Processing Toolbox algorithm ‘K-means clustering’ that can take a vector layer and group features into N clusters. A problem with this algorithm is that you do not have control over how many points end up in each cluster. Many applications require you to segment your data layer into equal sized clusters or clusters having a minimum number of points. Some examples where you may need this When planning for FTTH (Fiber-to-the-Home) network one may want to divide a neighborhood into clusters of at least 250 houses for placement of a node. Dividing a sales territory/ customers equally among sales teams with customers in the same region are assigned to the same team. There is a variation of the K-means algorithm called Constrained K-Means Clustering that uses graph theory to find optimal clusters with a user supplied minimum number of points belonging to given clusters. Stanislaw Adaszewski has a nice Python implementation of this algorithm that I have adapted to be used as a Processing Toolbox algorithm in QGIS. Warning! I have heard feedback from users that this algorithm doesn’t work on all types of point distributions and may get stuck while finding an optimal solution. I am looking into ways to improve the code and will appreciate if you had feedback. …

Google Earth Engine (GEE) is a powerful cloud-based system for analysing massive amounts of remote sensing data. One area where Google Earth Engine shines is the ability to calculate time series of values extracted from a deep stack of imagery. While GEE is great at crunching numbers, it has limited cartographic capabilities. That’s where QGIS comes in. Using the Google Earth Engine Plugin for QGIS and Python, you can combine the computing power of GEE with the cartographic capabilities of QGIS. In this post, I will show how to write PyQGIS code to programmatically fetch time-series data, and render a map template to create an animated maps like below. …

If you have collected GPS tracks, you know that the results can have varying accuracy. The track points collected along a route are not always on the road and can be jittery. …

When trying to automate your GIS workflows, one important step is the production of maps. Creating and exporting maps in QGIS is done via the Print Layout. One can automate creation of maps via the a rich Python API using the QgsLayout class. See our post Rendering Print Layouts from QGIS Models for a no-code solution to exporting print layouts. …

When you want to buffer features that are spread across a large area (such as global layers), there is no suitable projection that can give you accurate results. This is the classic case for needing Geodesic Buffers - where the distances are measured on an ellipsoid or spherical globe. This post explains the basics of geodesic vs. planar buffers well. …

Google recently announced Earth Engine Noncommercial Tiers for Google Earth Engine. This is a big change that affects all the non-commercial users of the GEE platform. Before the change, if you qualified for non-commercial use, you could use Earth Engine without any restrictions or fees. With the introduction of non-commercial tiers, you now have monthly limits on how much compute you can use for free. Once you exceed the allocated monthly quota, your account will enter a Restricted mode that will slow down the computations triggered through your account. See our updated Sign-up for Google Earth Engine page with guidance on how to select your quota tier. If you are a non-commercial user of Earth Engine, you now need to monitor and manage your quota usage to ensure you comply with these limits. This post outlines the concepts and tools you can use for quota monitoring. The post covers the following topics: Understanding Earth Engine Compute Unit (EECU) Monitoring Quota using Google Cloud Console Monitoring Quota using Python (using Google Colab) …

In this post, you will learn how to work with the Open Buildings 2.5D Temporal data and download it for many useful downstream applications, such as Visibility Analysis, Population Modeling, and 3D Visualization. Google has two important large-scale AI-derived open building datasets: Open Buildings V3 Polygons: This was released a few years ago and contains all buildings polygons detected from Google’s corpus of high-resolution imagery. You can read more about it in our post Mapping Building Density with Open Building Datasets. Open Buildings 2.5D Temporal V1: This is a newer dataset that aims to extract useful attributes for buildings such as year of construction and building height. Since this data is derived from open-source medium-resolution Sentinel-2 imagery, it has temporal coverage from 2016-2023. A deep learning model was trained to predict building heights from Sentinel-2 images, so we also get the height information each year. There is a new global 2.5D building heights dataset called GlobalBuildingAtlas. See my tutorial on extracting this data using open-source cloud-native technologies. We will cover a Google Earth Engine workflow to process this data to make it usable in a GIS environment and extract a high-resolution Digital Surface Model (DSM). We will also see how to combine the Open Buildings V3 polygon building footprints with the Open Buildings 2.5D Temporal V1 data to create and extract yearly polygon datasets containing building heights that can be used in a GIS environment. Open Buildings 2.5D Temporal Data Combined with Open Buildings V3 Polygons and Visualized in QGS The post is divided into the following sections. Part 1: Extracting Building Height Raster and High-Resolution DSM Part 2: Extracting Building Footprints with Heights …

When exporting large rasters from Google Earth Engine, it is recommended that you split your exports into several smaller tiles. In this post, I will share the best practices for creating tiled exports in your target projection that can be mosaicked together without any pixel gaps or overlaps. They key concept is the use of the crsTransform to ensure that each individual tile is on the same pixel grid. …

A temporally consistent global multi-class time-series classification dataset is critical to understand and quantify long-term changes. Till now, the choices were limited to lower resolution datasets such as MODIS Landcover (2000-present) at 500m resolution or ESA CCI (1992-present) at 300m resolution. We now have a new dataset GLC_FCS30D that provides a high-resolution landcover time-series derived from the Landsat archive (1984-2022) at 30m resolution with 35 classes. This is a very valuable dataset for studying landscape dynamics at high resolution and the first of its kind to be available made available in the public domain. The source dataset was released on Zenodo and can be downloaded as GeoTIFF files. This data is also available in the Google Earth Engine Community catalog and can be used within GEE directly. In this post, I want to share some technical details and scripts to help you analyze this data using Google Earth Engine. You will learn How to access and pre-process the GLC_FCS30D dataset. How to visualize and compare landcover changes between 1985-2022. How to calculate landcover statistics and export a CSV with areas of each class for the entire time series over multiple regions. Update: Seeour guide to access GLC_GCS30D using STAC/XArray without dependency on GEE. …

In this post, we will learn how to build a regression model in Google Earth Engine and use it to estimate total above ground biomass using openly available Earth Observation datasets. A new version of this tutorial using the Satellite Embedding dataset is available at Regression with Satellite Embedding Dataset NASA’s Global Ecosystem Dynamics Investigation (GEDI) mission collects LIDAR measurements along ground transects at 30m spatial resolution at 60m intervals. The GEDI waveform captures the vertical distribution of vegetation structure and is used to derive estimates of Aboveground Biomass Density (AGBD) at each sample. These sample estimates of AGBD are useful but since they are point measurements - you cannot directly use them to calculate total aboveground biomass for a region. We can use other satellite datasets to build a regression model using the GEDI samples to map and quantify the biomass in a region. Regression Workflow This article shows how we can build and run the entire workflow in Earth Engine - from pre-processing to building a regression model to running the predictions. You will also learn some advanced techniques and best practices such as: How to fuse datasets of different resolutions by using setDefaultProjection() and reduceResolution() to align them to a common grid. How to sample pixels from rasters with sparse data efficiently and precisely by leveraging the image mask using stratifiedSample(). How to split your workflow into separate steps and use Exports to avoid user memory limit exceeded or computation timed-out errors. …

Dynamic World is a new landcover product developed by Google and World Resources Institute (WRI). It is a unique dataset that is designed to make it easy for users to develop locally relevant landcover classification easily. Contrary to other landcover products which try to classify the pixels into a single class - the Dynamic World (DW) model gives you the probability of the pixel belonging to each of the 9 different landcover classes. The full dataset contains the DW class probabilities for every Sentinel-2 scene since 2015 having <35% cloud-cover. It is also updated continuously with detections from new Sentinel-2 scenes as soon as they are available. This makes DW ideal for change detection and monitoring applications. A key fact about this dataset is that Dynamic World is not a ready-to-use landcover product. Users are expected to fine-tune the output of DW with local knowledge into a final landcover product. Since DW provides per-pixel probabilities generated by a Fully Convolutional Neural Network (FCNN) model, a lot of difficult problems encountered in classifying remotely sensed imagery are addressed already and allows users to refine it with a relatively simple model (such as Random Forest) with small amount of local training data. A good mental model to use for Dynamic World is to not think of it as landcover product but as a dataset that provides 9 additional bands of landcover related information for each Sentinel-2 image that can be refined to build a locally relevant classification or change detection model. As seen in the mangrove classification example, using the Dynamic World probability bands as input to a supervised classification model can help you generate a more accurate landcover map in less amount of time. It also eliminates the need for post-processing the results. To test this concept and explore the potential of this new dataset in developing locally relevant landcover maps - I partnered with Google and WRI to develop a training workshop and host a 5-day “Mapathon” with participants of diverse backgrounds. The event was a mix of hands-on workshop along with hackathon-style group projects to use Dynamic World for a real-world application. The workshop was hosted by Regional Centre for Mapping of Resources for Development (RCMRD) in Nairobi, Kenya. You can read more about the event in this article. I and Elise Mazur from WRI also gave a talk about our experience at Geo for Good 2023. In this post, I want to share more technical details about the workshop materials and code for projects for those who may want to use Dynamic World for their own applications. …

Extracting building footprints from high-resolution imagery is a challenging task. Fortunately we now have access to ready-to-use building footprints dataset extracted using state-of-the-art ML techniques. Google’s Open Buildings project has mapped and extracted 1.8 billion buildings in Africa, South Asia, South-East Asia, Latin America and the Caribbean. Microsoft’s Global ML Building Footprints project has made over 1.24 billion building footprints available from most regions of the world. Update: VIDA has released the most comprehensive buildings dataset by combining both Google and Microsoft building footprint dataset. This dataset is available in Google Earth Engine via the GEE Community Catalog. Update: Google has released another dataset for building heights. See our post on Extracting Building Heights from Open Buildings 2.5D Temporal Dataset Given the availability of these datasets, we can now analyze them to create derivative products. In this post, we will learn how to access these datasets and compute the aggregate count of buildings within a regular grid using Google Earth Engine. We will then export the grid as a shapefile and create a building density map in QGIS. …

An important concept in spatial statistics is pixel weights. When calculating pixel statistics with a polygon, partial pixel overlaps are treated differently by different packages and you need to understand this to evaluate the accuracy of your results. Consider the following image. What is the correct answer? …

In this article, I will outline a method for extracting shoreline from satellite images in Google Earth Engine. This method is scalable and automatically extracts the coastline as a vector polyline. The full code link is available at the end of the post. UPDATE: The post now includes tidal-phase filtering using HYCOM Data. The method involves the following steps Create a Cloud-free Composite Image from images collected during the same tidal phase Extract All Waterbodies Remove Inland Water and Small Islands Convert Raster to Vector Simplify and Extract Coastline Video Demonstration of the Script We will go through the details of each step and review the Google Earth Engine API code required to achieve the results. …

If you are like me, you have a lot of assets uploaded to Earth Engine. As you upload more and more assets, managing this data becomes quite a cumbersome task. Earth Engine provides a handy Command-Line Tool that helps with asset management. While the command-line tool is very useful, it falls short when it comes to bulk data management tasks. What if you want to rename an ImageCollection? You will need to manually move each child image to a new collection. If you wanted to delete assets matching certain keywords, you’ll need to write a custom shell script. If you are running low on your asset quota and want to delete large assets, there is no direct way to list large assets. Fortunately, the Earth Engine Python Client API comes with a handy ee.data module that we can leverage to write custom scripts. In this post, I will cover the following use cases with full python scripts that can be used by anyone to manage their assets: How to get a list of all your assets (including folders/sub-folders/collections) How to share all assets in a folder How to find the quota consumed by each asset and find large assets How to rename ImageCollections How to delete ImageCollections The post explains each use-case with code snippets. If you want to just grab the scripts, they are linked at the end of the post. …

Many applications require replacing missing pixels in an image with an interpolated value from its temporal neighbours. This gap-filling technique is used in several applications, including: Replacing Cloudy Pixels: You may want to fill gap in an image with the best-estimated value from the before and after cloud-free pixel. Estimating Intermediate Values: You can use this technique to compute an image for a previously unknown time-step. If you had population rasters at 2 different years and want to compute a population raster for an intermediate year using pixel-wise linear interpolation. Preparing Data for Regression: All of your independent variables may not be available at the same temporal resolution. You can harmonize various dataset by generating interpolated raster at uniform or fixed time-steps. Google Earth Engine can be used effectively for gap-filling time-series datasets. While the logic for linear interpolation is fairly straightforward, data preparation for this in GEE can be quite challenging. It involves use of Joins, Masks and Advanced Filters. This post explains the steps with code snippets and builds a fully functioning script that can be applied on any time-series data. You can also watch my video at Geo for Good 2022 Conference where I covered the same topic. …

Most optical satellite imagery products come with one or more QA-bands that allows the user to assess quality of each pixel and extract pixels that meet their requirements. The most common application for QA-bands is to extract information about cloudy pixels and mask them. But the QA bands contain a wealth of other information that can help you remove low quality data from your analysis. Typically the information contained in QA bands is stored as Bitwise Flags. In this post, I will cover basic concepts related to Bitwise operations and how to extract and mask with specific quality indicators using Bitmasks. …

In this post, I will outline techniques for computing weighted-centroids in both QGIS and Google Earth Engine. For a polygon feature, the centroid is the geometric center. It can also be thought of as the average coordinate of all points within the polygon. There are some uses cases where you may want to compute a weighted-centroid where some parts of the polygon gets higher ‘weight’ than others. The main use-case is to calculate a population-weighted centroid. One can also use Night Lights data as a proxy for urbanized population and calculate a nightlights-weighted centroid. Some applications include: Regional Planning: Locate the population-weighted centroid to know the most accessible location from the region. Network Analysis: For generating demand points in location-allocation analysis, you need to convert demands from regions to points. It preferable to compute populated-weighted centroids for a more accurate analysis. Do check out this twitter-thread by Raj Bhagat P for more discussion on weighted centroids. Different Weighted Centroids for the State of Karnataka, India (2015) …

Google Earth Engine makes it easy to compute statistics on gridded raster datasets. While calculating statistics on imagery datasets is easy, special care must be taken when working with population datasets. In this post, I will outline the correct technique for computing statistics for population rasters and aggregating pixels. …

Many useful climate and weather datasets come as gridded rasters. The techniques for working with them is slightly different than other remote sensing datasets. In this post, I will show how to work with gridded rainfall data in Google Earth Engine. This post also serves an an example of how to use the map/reduce programming style to efficiently work with such large datasets. …

Color correction is an important process working with satellite and aerial imagery. A common technique used to balance the colors across multiple images is Histogram Matching. While the algorithm has been around for a long time, there aren’t many free and open-source tools that can used at scale. Mapbox has released an open-source tool called rio-hist that works well for small and medium sized images. Whitebox Tools has a Histogram Matching algorithm that can be used in QGIS via Whitebox Tools Processing Plugin. But when working with large mosaics, such as the ones used in this post - it runs out of memory or takes a very long time. Google Earth Engine is a good alternative to perform fast histogram matching across large images. In this post, I will first give an overview of the histogram matching algorithm and then show you how it can be implemented in Earth Engine. The example images are large high resolution orthomosaics (3cm/pixel resolution) collected by UAV around Oakland, CA area. Source Images: © Dan Koopman, Downloaded from OpenAerialMap …

When working on Remote Sensing applications, many operations require calculating area. For example, one needs to calculate area covered by each class after supervised classification or find out how much area within a region is affected after a disaster. Calculating area for rasters and vectors is a straightforward operation in most software packages, but it is done in a slightly different way in Google Earth Engine - which can be confusing to beginners. In this post I will outline methods of calculating areas for both vectors as well as images. We will cover the following topics, starting from simple to complex. Area Calculation for Features i.e. vector data Area Calculation for Images (Single Class) Area Calculation for Images by Class Area Calculation for Images by Class by Region Area Calculation for Images by Class by Region by Year …

Time series analysis is one of the most common operations in Remote Sensing. It helps understanding and modeling of seasonal patterns as well as monitoring of land cover changes. Earth Engine is uniquely suited to allow extraction of dense time series over long periods of time. In this post, I will go through different methods and approaches for time series extraction. While there are plenty of examples available that show how to extract a time series for a single location - there are unique challenges that come up when you need a time series for many locations spanning a large area. I will explain those challenges and present code samples to solve them. The ultimate goal for this exercise is to extract NDVI time series from Sentinel-2 data over 1 year for 100 farm locations spanning an entire state in India. A new version of this tutorial using open-source cloud-native technologies is available at Extracting NDVI Time-Series for Multiple Field Boundaries. This tutorial shows how to accomplish the same without dependency on GEE. …

Google Earth Engine (GEE) is a powerful cloud-based system for analysing massive amounts of remote sensing data. One area where Google Earth Engine shines is the ability to calculate time series of values extracted from a deep stack of imagery. While GEE is great at crunching numbers, it has limited cartographic capabilities. That’s where QGIS comes in. Using the Google Earth Engine Plugin for QGIS and Python, you can combine the computing power of GEE with the cartographic capabilities of QGIS. In this post, I will show how to write PyQGIS code to programmatically fetch time-series data, and render a map template to create an animated maps like below. …

Historically, the open-source spatial analysis workflows lived in different silos Python: You use Pandas, GeoPandas, and Jupyter Notebooks. SQL: You use SQL with PostGIS. If you lived in the Python ecosystem, you rarely have to switch to SQL for analysis and vice versa. Many people, including me, had little motivation to sharpen their SQL skills when you could get away with doing things in Python. Many of our students who wanted to learn SQL, found themselves choosing between these two stacks and found SQL had a much higher friction to get started. Recently, I have been using DuckDB and find that it is the perfect bridge between these two ecosystems. DuckDB + Python + LLMs provide the easiest and most rewarding pathway for Python users to learn and incorporate SQL in their workflows. In this post, we will cover the following topics What is DuckDB? Using DuckDB to learn SQL with the help of LLMs Example Workflow with LLMs Querying and Loading Administrative Boundaries from GeoBoundaries Advanced Workflows Extracting Overture Maps Data Extracting Farm Boundaries from Global Fields of The World (FTW) Open our companion notebook learning_sql_duckdb.ipynb in Google Colab to follow along and run the queries yourself. …

Welcome to #PythonDatavizChallenge - Learn Mapping and Data Visualization with Python in 30 Days! We have designed this challenge to help you learn how to create charts, maps, animations, dashboards and interactive mapping applications using Python ! Spend 30 minutes each day for the next 30 days to level-up your Python dataviz skills. We have spent over 2 years building and refining this course and are excited to share it with you all - completely free. The challenge is a series of short videos, one set for each day, that cover the full course material step by step. The material covers both static and dynamic plotting libraries along with the app framework - Streamlit. At the end of the course, you will have the necessary skills to build data-powered web mapping apps and dashboards. Ready for #PythonDatavizChallenge? Read on to know the details. This is an intermediate course that assumes good working knowledge of Python. If you are new to programming, complete our Python Foundation for Spatial Analysis course first. …

ISRO recently released the full archive of medium and low-resolution Earth Observation dataset to the public. This includes the imagery from LISS-IV camera aboard ResourceSat-2 and ResourceSat-2A satellites. This is currently the highest spatial resolution imagery available in the public domain for India. In this post, I want to cover the steps required to download the imagery and apply the pre-processing steps required to make this data ready for analysis - specifically how to programmatically convert the DN values to TOA Reflectance. We will use modern Python libraries such as XArray, rioxarray, and dask - which allow use to seamlessly work with large datasets and use all the available compute power on your machine. …

An important concept in spatial statistics is pixel weights. When calculating pixel statistics with a polygon, partial pixel overlaps are treated differently by different packages and you need to understand this to evaluate the accuracy of your results. Consider the following image. What is the correct answer? …

Matplotlib has functionality to created animations and can be used to create dynamic visualizations. In this post, I will explain the concepts and techniques for creating animated charts using Python and Matplotlib. I find this technique very helpful in creating animations showing how certain algorithms work. This post also contains Python implementations of two common geometry simplification algorithms and they will used to create animations showing each step of the algorithm. Since both of these implementations use a recursive function, the technique shown in the post can be extended to visualize other recursive functions using matplotlib. You will learn how to create animated plots like below. …

Spatial indexing methods help speed up spatial queries. Most GIS software and databases provide a mechanism to compute and use spatial index for your data layers. QGIS as well as PostGIS use a spatial indexing scheme based on R-Tree data structure - which creates a hierarchical tree using bounding boxes of geometries. This is quite efficient and results in big speedup in certain types of spatial queries. Check out Spatial Indexing section of my course Advanced QGIS where I show how to use R-Tree based Spatial index in QGIS. If you use Python for geoprocesisng, the GeoPandas library also provides an easy to use implementation of R-Tree based spatial index using the .sidex attribute. University of Helsinki’s AutoGIS course has an excellent example of using spatial index with geopandas. In this post, I want to talk about another spatial indexing system called H3. …

When working with raster data, you may sometimes need to deal with data gaps. These could be the result of sensor malfunction, processing errors or data corruption. Below is an example of data gap (i.e. no data values) in aerial imagery. Source Image: © Commission for Lands (COLA) ; Revolutionary Government of Zanzibar (RGoZ), Downloaded from OpenAerialMap. (Note: The data gap is simulated using a python script and is not part of the original dataset) …

When working with raster datasets of different projections and resolutions, it is often desirable to reproject them to the same projection and align them to the same pixel grid. In this post, we will explore the recently introduced options in the open-source GDAL utility gdalwarp that makes this process much simpler and efficient. In particular, we will be exploring the -r sum (Resample with Sum), -r average (Resample with Average) and -tap (Target Aligned Pixels). We will take the following 3 raster datasets and clip, resample and align them to a common pixel grid. LandScan Global: A high-quality global population grid that is available at 1km resolution in the geographic CRS WGS84 Lat/Lon (EPSG:4326). GHS Population Grid: A 100m resolution global population dataset that is distributed in the World Mollweide Equal Area Projection (ESRI:54009). NLCD Tree Canopy Cover: A 30m resolution gridded dataset with percent canopy estimate of tree cover in the NAD83 CONUS Alberts Projection (EPSG:5070). As you can see we have datasets that have widely varying pixel sizes and projections. If we wanted to compare them with each other - we must first harmonize them on a unified pixel grid. We will learn how to reproject, resample and align these to the NAD83 California Albers Projection (EPSG:3311) and at 1km resolution. …

GDAL has support for GDAL Virtual File System which allows GDAL library and command-line tools to work with files on network storage. This is critical in the era of Cloud-Native Geospatialwhere it is becoming a standard practice to access and share geospatial data via cloud storage services. In this post, we will see how to use GDAL command-line tools to read and write data to Google Cloud Storage (GCS) using the /vsigs file system handler. We will focus exclusively on Google Cloud Storage for this post -but the same concepts apply when using other cloud services such as AWS S3 or Azure Blobstore. Similarly, other GDAL-based tools - such as rasterio - will be able to access the data from GCS using the same configuration options shown here. The post covers the following topics Reading Files from Public GCS Buckets. Creating Private GCS Buckets and Uploading Data Configuring Authentication and Reading Data from Private GCS Buckets Writing Data to Private GCS Buckets Using Environment Variables This post assumes familiarity with the GDAL command-line tools and assumes you have installed GDAL on your machine. You will find detailed instructions for installation in our course material for Mastering GDAL Tools. The code snippets are split over multiple lines for redability using the Windows Line Continuation character ^. If you are running these on Mac/Linux, you can replace them with \ instead. …

Multi-criteria Overlay Analysis is the process of the allocation of land to suit a specific objective on the basis of a variety of attributes that the selected areas should possess. Although this is a common GIS operation, it is best performed in the raster space. This post outlines the typical workflow to take source vector data, transform them to appropriate rasters, re-classify them and perform mathematical operations to do a weighted suitability analysis. The post uses the command-line utilities provided by the open-source GDAL library. If you want to do such analysis in QGIS, please check my tutorial at Multi Criteria Overlay Analysis (QGIS3) We will work with crime and infrastructure data for the city of London and find suitable areas to build new parking facilities that can help reduce bicycle thefts. Our analysis will apply the following 3 criteria. The proposed parking must be In a bicycle theft hotspot Close to a bicycle route Far from existing parking facilities The problem statement: Identify suitable locations for building new bicycle parking facilities …

When working with raster data, you may sometimes need to deal with data gaps. These could be the result of sensor malfunction, processing errors or data corruption. Below is an example of data gap (i.e. no data values) in aerial imagery. Source Image: © Commission for Lands (COLA) ; Revolutionary Government of Zanzibar (RGoZ), Downloaded from OpenAerialMap. (Note: The data gap is simulated using a python script and is not part of the original dataset) …

gdal_calc is one of the lesser used tools among the GDAL utilities. There aren’t many examples of using it in the wild and some advanced features are not well documented. I am finding myself using it a lot lately and have discovered some really powerful use cases. …

When working with elevation data, sometimes you may discover that 2 datasets from different providers have very different elevation values for the same location. A common reason for this being each dataset being referenced to a different surface. …

GeoPDF is a unique data format that brings the portability of PDF to geospatial data. A GeoPDF document can present raster and vector data and preserve the georeference information. This can be a useful format for non-GIS folks to consume GIS data without needing GIS-software. While GeoPDF is a proprietary format, we have a close alternative in the open Geospatial PDF format. GDAL has added support for creating Geospatial PDF documents from version 1.10 onwards. In this post, I will show how to create a GeoPDF document containing multiple vector layers. …

GDAL and OGR libraries come with handy command-line tools. These tools are quite powerful and can save you a lot of effort if you know how to use them. Here I will show how to use the ogrinfo and ogr2ogr tools to perform spatial joins. A single command can do complex operations on your spatial data and save you a lot of clicking-around and data-munging in a GIS. …

In a previous post, I showed how to use the Uber Movement Travel Times data to create isochrones. In this post, we will explore another use case of this dataset. Say you are concerned about loss of productivity due to long commute times of your employees and wonder if a change in office times might help them get to the office faster. A similar analysis can be done to see if a change in office location will result in better or worse commutes. Here’s the hypothetical scenario “Given the location of an office and location of homes of employees, determine their current commute times for office timings of 9am-11am and 5pm-6pm. If the office timings were changed to non-peak timings of 7am-8am and 3pm-4pm, what would be the time savings?” …

Uber Movement data has been discontinued. This post now serves as a reference only for similar analysis. Mapshaper is a free and open-source software for spatial data processing. It is written in javascript and runs in your browser without any extra plugins and can perform a range of analysis. It started out as a tool for topologically-aware simplification, but has evolved into a swiss army knife of spatial data processing tools. All processing is done in the browser locally, and I have found that it can handle large volumes of data easily and processing is usually much faster than desktop based GIS software. …

Mapshaper is a free and open-source tool that is best known for fast and easy simplification. Other tools for simplification - like QGIS or ogr2ogr - do not preserve topology while simplifying. This means you may get sliver polygons or missing intersections. Mapshaper performs topologically-aware simplification and gives you much more control on the process. …

We recently moved spatialthoughts.com from WordPress to a static site built with Hugo. The new site is a static site that is deployed with GitHub Actions to GitHub Pages. The redesign is not merely a frontend change, but enables us to automate the whole process of running our live cohort-based classes. The interesting part of the whole migration is the backend design - where all the data for our courses and cohorts is stored as YAML files - and all pages are generated from them. This post explains our design choices and how the new site is put together from the data files. The new spatialthoughts.com home page, built with Hugo …

I was recently asked to deliver a session on Earth Science for kids. My daughter goes to an after-school science program at Max Science where they teach science with unique and fun hands-on experiments. They wanted to do an interactive session to introduce the kids to Earth Science and asked if I could deliver a guest talk for their kids in Grades 1 to 4. I loved the idea and developed a module titled The Science of Satellites to introduce the magic of remote sensing to primary school kids. The session ended up being a lot of fun, for the kids and me. In this post, I want to go through the materials and my experience teaching this session. All the content developed for this session - including high-resolution graphics - is available freely for download. Scroll to the bottom to find the download link. The 1.5 hour session was split into 3 parts: Part 1 Guess the Place A game to guess the place from satellite images Part 2: The Science of Satellites Learning what satellites do and how can you build and launch a satellite Part 3: Your Name from Space An activity where kids create their from letters seen from satellite images. …

As everyone who is involved in teaching and training knows, the past few months have been hard. We all had to make changes to accommodate working from home and adopting online teaching methods. Before the COVID-19 outbreak, I used to conduct all my training in-person. Either hosting it at a training center or at a client location. My materials, structure and instruction style was tuned to this setup. I was skeptical whether the experience of a classroom can be replicated - even partially - online. Over the past 2 months, I have conducted numerous online training sessions. All my courses have been moved to a ’live’ online class and even started offering short-format classes. I did a lot of research, talked to other trainers and spent a considerable effort in trying to make this transition. I thought sharing some of the lessons and best practices here will help fellow educators. …

I was invited to participate in a panel discussion on Geospatial Intelligence for #LetsTalkDeepTech Webcast hosted by Swiggy. I talked about the history and evolution of this space and gave a deep dive into solutions for deriving intelligence from imagery. Below is the a longer version of my talk on evolution of location intelligence with some references. I also share a copy of my presentation at the end. Hope you find it useful. Agree/Disagree with my views? Let me know in the comments. …

Converting between different coordinate reference systems or projection is a fairly standard feature in a GIS. There are a number of command line tools also available for performing bulk-conversions. cs2cs program, part of the PROJ.4 library is my favorite. …

This weekend, I got an opportunity to volunteer with a non-profit called Junglescapes. We took a day trip to the Bandipur forest in Karnataka where they have done extensive work in forest restoration. One of their success stories is working with the locals to remove invasive species such as Lantana from the forest. Junglescapes volunteers and locals carry out regular line transact surveys to determine the impact of their interventions. One of the goals for my participation was to see if we can replace the cumbersone paper forms and handheld GPS devices with a mobile-phone based survey using ODK. I am sharing my notes on how we setup the survey and mapping of the result. …

The Association for People with Disability (APD) is a non-profit organization based out of Bangalore, India. Their mission is to reach out and rehabilitate people with disability from the under privileged segment. Over the past year, I along with my colleagues have been volunteering with them to develop a system that can help improve their field data collection efforts. …

A long pending weekend project is done. Printed, cut and folded a sturdy globe using the template from Le Paper Globe. This is not only fun, but a good prop to learn more about Geography. I envision it would make a fun do-it-yourself project with kids of all ages. …

I recently had a need to calculate distance between a large number of latitude/longitude coordinate pairs. There are many options available if you want to import these in a GIS and run analysis. But there is a simpler and much more accesible way if you aren’t doing very high accuracy calculations. …

Many students have asked me for ideas on what topic they should choose for their Thesis. I have debated this myself when I was a student. The ideal topic would be the one that allows you to dive into a topic deeply as well as give you some practical skills that will help you landing a job. Here are some pointer that may help you make that decision. …