محاسبه منطقه ای میانگین بارش(QGIS3)¶

محاسبه بیلان آب، مدل سازی سیلاب، پیش بینی رواناب، مطالعات اقلیمی و غیره اغلب به عمق متوسط بارندگی در یک حوضه هیدرولوژیکی به عنوان ورودی نیاز دارد - که همچنین به نام بارش منطقه ای یا میانگین بارندگی (AMR) نامیده می شود.

محاسبه AMR را می توان با استفاده از داده های باران سنج انجام داد. با استفاده از محل بارش باران و بارش مشاهده شده، می توان با استفاده از هر یک از تکنیک های زیر، بارش متوسط را در یک مکان مشخص تخمین زد:

** میانگین ریاضی **: می توان به سادگی به طور متوسط از تمام مقادیر مشاهده شده استفاده کرد. این روش فرض می کند که میدان بارندگی همگن است و مشاهدات در باران نگارها مستقل است و وزن برابری به تمام سنجه های باران می دهد.
** چندضلعی تیسن Thiessen **: این روش منطقه را با استفاده از چند ضلعی Thiessen تقسیم می کند با این فرض که بارندگی در محدوده پوشش هر ایستگاه همگن است. این روش همچنین به طور متوسط میانگین وزن منطقه ای نامیده می شود. این مفروضات برای زمین های کم ارتفاع و یا مسطح خوب است، اما برای زمین های کوهستانی مناسب نیست.
** روش Iso-hyetal *: این تکنیک درونیابی خطوط هم تراز Isohyets - خطوط با بارش یکسان با محاسبه می کند. فرض بر این است که بارندگی بین دو خط تراز Isohyets همگن است. این روش زمانی مناسب است که شبکه باران سنجی متراکم است.
** فاصله وزنیWeighting / شبکه Gridded ** - این یک تکنیک درون یابی است که در آن یک شبکه شطرنجی ایجاد می شود و مقدار برای هر پیکسل بر اساس فاصله تا ایستگاه تخمین زده می شود. هنگامی که نقاط شبکه تخمین زده شده اند، جمع شده اند و مجموع با تعداد سلول های شبکه تقسیم می شود تا میانگین بارش منطقه ای را بدست آورد.
** روش های زمین آماری Geostatistic**: بارندگی به شدت تحت تأثیر عوامل محلی - مانند ارتفاع قرار دارد. با استفاده از تکنیک های رگرسیون چند متغیره یا تکنیک های کریگینگ، می توان ازخودهمبستگی مکانی حساب کرد و می تواند دقت بیشتری به دست آورد. این روش ها زمانی مناسب هستند که توزیع ایستگاه سنجنده باران یکنواخت و متراکم باشد.

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.

توجه

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.