Calculating Areal Mean Rainfall (QGIS3)¶

Calculation of water balance, flood modeling, runoff forecasting, climate studies etc. often need the average depth of rainfall in a hydrological basin as an input - which is also called Areal Precipitation or Areal Mean Rainfall (AMR).

AMR calculation can be done using rain gauge data. By using the rain gauge location and observed precipitation, one can estimate the average precipitation at a given location by using any of the following techniques:

Arithmetic Average: One can simply take an average of all the observed values. This method assumes that the rainfall field is homogeneous and that the rain gauge observations are independent and give equal weight to all rain gauges.
Thiessen Polygon: This method divides the area using Thiessen polygons with the assumption that rainfall is homogeneous within the coverage area of each station. This method is also called an area-weighted average. These assumptions are fine for low-lying or flat terrain, but not suitable for mountainous terrain.
Iso-hyetal Method: This interpolation technique calculates Isohyets - lines joining equal precipitation. It assumes that rainfall between 2 isohyets is homogeneous. This method is suitable when the rain-gauge network is dense.
Distance Weighting/Gridded - This is an interpolation technique where a raster grid is created and a value for each pixel is estimated based on the distance to stations. Once the grid points have all been estimated they are summed and the sum is divided by the number of grid cells to obtain the areal mean precipitation.
Geostatistical Methods: Rainfall is strongly influenced by local factors - such as elevation. Using multivariate regression or Kriging techniques, one can account for spatial autocorrelation and can achieve better accuracy. These methods are suited when the distribution of the rain gauge station is uniform and dense.

In this tutorial, we will learn the QGIS workflow to calculate Areal Mean Rainfall using the Thiessen Polygon method.

Overview of the task¶

In this tutorial, we will take the precipitation measured by the Global Historical Climatology Network (GHCN) stations and compute Areal Mean Rainfall in each Hydrological basin in the state of Florida.

Other skills you will learn¶

How to remove data with Null values.
How to fix invalid geometries in a layer.
How to check your Processing History and re-run a tool with the same parameters.
How to dissolve polygons and summarize statistics.
How to use only selected features in Processing algorithms.

Get the data¶

We will use NOAA Climate data , HydroSHEDS and US Census Bureau Cartographic Boundary data layers.

Station-wise Precipitation¶

Go to NOAA Climate data website. Click on the Search tool.

Select Global Summary of the Month in Select weather Observation Type/Dataset, then in Select Date Range choose July 2020, in Search For select States, in Enter a Search Term type Florida. Click Search.

Click Add To Cart, in items click View Items.

Switch to Custom Global Summary of the Month CSV, click Continue.

Check the Geographic Locations and in Select data types for custom output expand Precipitation select Precipitation (PRCP). Click Continue.

Enter the mail address and click SUBMIT ORDER to Download the data.

Hydrological Basins¶

Go-to HydroSHEDS website and click Download.

Select HydroBASINS ‣ Standard ‣ North America and Caribbean ‣ hybas_na_lev06_v1c.zip

Enter the mail address and click Submit Request to Download the data.

State Boundaries¶

Visit the Cartographic Boundary Files - Shapefile page. Download the cb_2018_us_state_500k.zip file from the States section.

For convenience, you may directly download a copy of all the datasets from the links below:

florida_2020_07_prcp

hybas_na_lev06_v1c

cb_2018_us_state_500k

Data Sources: [GHCN], [HYDROSHEDS], [USCENSUS]

Procedure¶

Open QGIS and click on the Open Data Source Manager.

In the Data Source Manager dialog box, switch to Delimited Text. Click on the … in File name then browse and select the florida_2020_07_prcp.csv file.

Now, under Geometry Definition choose Point coordinates, X field and Y field should be Longitude and Latitude respectively. Choose the Geometry CRS as EPSG 4326 - WGS 84. Click Add.

Now a new point layer will be added, click on the Open Attribute Table icon.

In the Attribute table the field PRCP represents the amount of precipitation in the station during the July 2020. Also, this data is recorded in inches. Note there are few Null values which can cause problems during calculation. Sort the PRCP column, and you would see there is only a small fraction of the dataset is Null. We will now remove the stations with Null values.

Open the Processing Toolbox by going to Processing ‣ Toolbox, and search and locate the Vector selection ‣ Extract by attribute algorithm.

In the Extract by Attribute dialog box, Select the Input layer as florida_2020_07_prcp, then choose PRCP in Selection attribute, then is not null in Operator. Click on the … next to Extracted (attribute), choose Save to File..., enter the layer name as precipitation_filtered.gpkg and click Run.

Now a new layer is added to canvas, turn off the older layer, and you can see the stations with Null values have been removed.

Now we will generate thiessen polygons from this layer. Open the processing toolbox by going to Processing ‣ Toolbox, and search and locate the Vector Geometry ‣ Voronoi polygon algorithm.

Catatan

Thiessen polygons represent the area of influence around each point. Every polygon defines the area which is closer to a particular station than any other station. This means the precipitation at any point is assumed to be the same as the nearest station.

Select precipitation_filtered as the Input layer. Since we do not have data for any rain-gauge stations outside the basin boundary, we can add some buffer area. Enter 15 as the Buffer region (% of extent). Click on the … in Voronoi polygons and select Save to File…, then enter the name as thiessen_polygons.gpkg. Click Run.

A new layer thiessen_polygons will be added to canvas. Let's clip this layer to the state boundary. Search for cb_2018_us_state_500k.shp file in Browser and drag it to canvas.

The states layer is in a different CRS than the Project CRS. You will get a prompt with different options for transforming this CRS to the Project CRS. In Select Transformation Dialog box, you can choose the default transformation and click OK.

We will now clip the thiessen polygons layer to the Florida state boundary. Click on the Select Feature by area or Single Click icon and click over Florida state.

Open the Processing Toolbox by going to Processing ‣ Toolbox, and search and locate the Vector overlay ‣ Clip algorithm.

In the Clip dialog box, select the Input layer as thiessen_polygons, in the Overlay layer select the cb_2018_us_state_500K layer and check the Selected features only checkbox, then click on the … in Clipped and select Save to File… , then enter the name as thiessen_polygons_clipped.gpkg. Click Run.

The clipped thiessen polygons layer thiessen_polygons_clipped will be added to the canvas. Turn-off the visibility of all other layers. As our task is to calculate average rainfall over each basin, we will now load the polygons representing basins. Locate the hybas_na_lev06_v1c.shp layer from the Browser and add it to the canvas.

You will notice that each basin is covered by many thiessen polygons and each polygon spans multiple basins. To visualise this Open layer styling panel icon and change the Opacity to 75%. We will now intersect both the layers to cut the thiessen polygons to the boundary of each basin.

Open the Processing Toolbox by going to Processing ‣ Toolbox, and search and locate the Vector overlay ‣ Intersection algorithm.

In the Intersection dialog box, select the Input layer as thiessen_polygons_clipped and Overlay layer as hybas_na_lev06_v1c, then click on the … in Intersected and select Save to File… , then enter the name as thiessen_polygons_basin.gpkg. Click Run.

The execution will fail with an error message has invalid geometry. Please fix the geometry or change the Processing setting to the “Ignore invalid input features” option.. You can learn more about this error in the Handling Invalid Geometries (QGIS3) tutorial.

To fix the geometries, open the Processing Toolbox by going to Processing ‣ Toolbox, and search and locate the Vector geometry ‣ Fix geometries algorithm.

In the Fix Geometries dialog box select the Input layer as hybas_na_lev06_v1c and click on … on Fixed geometries and select Save to File…, enter the file name as hybas_na_lev06_v1c_fixed.gpkg and click Run.

Now a new layer will be added to canvas. We can now try the intersection again. Instead of running the tool from scratch and filling all the parameters, we can retrieve the pre-filled dialog from Processing History and modify only the Overlay layer. Click Processing ‣ History.

Double-click on the native:intersection algorithm from the list.

Change the Overlay layer to hybas_na_lev06_v1c_fixed and click Run.

Now a new layer will be loaded, and you can see the thiessen_polygons_basin is clipped based on the basin boundary.

Now, let’s calculate the average rainfall value from the thiessen polygons for each basin. This is done using the Aggregate tool which allows us to dissolve individual polygons while calculating statistics on the attribute values. Now, open the Processing Toolbox by going to Processing ‣ Toolbox, and search and locate the Vector geometry ‣ Aggregate algorithm.

In the Aggregate dialog box choose Input layer as thiessen_polygons_basin, select all fields except PRCP and HYBAS_ID and click Delete selected field.

In Group by expression select HYBAS_ID. This means that the tool will dissolve all polygons that have the same value of HYBAS_ID. In our case, these will be all thiessen polygons falling a basin. In the Aggregates section, we can configure how different field values will be aggregated from all polygons that gets dissolved. For PRCP, click on the expression button to enter the below expression. The expression calculates the area-weighted fraction for each polygon. Set the Aggregate Function to sum, which will sum up all the area-weighted fractions resulting in the area-weighted mean. For HYBAS_ID, change the Aggregate Function to first_value. Since we are grouping all thiessen polygons by their HYBAS_ID, all the values will be the same and the first_value function will use the attribute value from the first polygon in each basin. Click on … on Aggregated and select the Save to File…, enter the file name as areal_mean_rainfall.gpkg and click Run.

(PRCP * $area) / sum($area)

A new layer will be added to canvas, lets open the Attribute table to explore. Click on the Open Attribute Table icon.

The field PRCP contains the areal mean rainfall for each basin in inches.

If you want to report any issues with this tutorial, please comment below. (requires GitHub account)