# 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:

1. 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.

2. 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.

3. 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.

4. 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.

5. 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.

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.

## Alte competențe pe care le veți dobândi¶

• 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.

## Obținerea datelor¶

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

### Station-wise Precipitation¶

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

1. 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.

1. Click Add To Cart, in items click View Items.

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

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

### Hydrological Basins¶

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

### State Boundaries¶

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

Pentru comoditate, puteți descărca o copie a seturilor de date direct de la adresele de mai jos:

florida_2020_07_prcp

hybas_na_lev06_v1c

cb_2018_us_state_500k

Sursele datelor: [GHCN], [HYDROSHEDS], [USCENSUS]

## Procedura¶

1. Open QGIS and click on the Open Data Source Manager.

1. 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.

1. 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.

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

1. 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.

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

1. 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.

1. 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.

1. 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.

Notă

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.

1. 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.

1. 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.

1. 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.

1. 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.

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

1. 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.

1. 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.

1. 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.

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

1. 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.

1. 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 Gestionarea Geometriilor Nevalide (QGIS3) tutorial.

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

1. 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.

1. 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.

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

1. Change the Overlay layer to `hybas_na_lev06_v1c_fixed` and click Run.

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

1. 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.

1. 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.

1. 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)
```
1. A new layer will be added to canvas, lets open the Attribute table to explore. Click on the Open Attribute Table icon.

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

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)