Repository navigation
How to generalise specifying spatial information? #66
Description
So far, HydPy misses a standard way to define spatial information. Here are some first ideas.
Where do we need spatial information so far?
L-Land Version 3 requires the latitutde and the longitude in decimal degrees for calculating extraterrestrial radiation. Each spatial location corresponds to the area centre of a sub-catchment. Each subcatchment is represented by at least one Element object.
The same holds for Evap Version 1.
All Conv models need information on multiple locations. Parameter InputCoordinates defines the coordinates of all relevant input points, for example of different rain gauges. Parameter OutputCoordinates defines the coordinates of all relevant output points, for example of the area centres of all relevant sub-catchments. The unit is arbitrary, but accurate interpolations usually require a projection system like UTM. Each input or output location is mapped to a Node object.
Conv Version 3 also requires additional information for each input and output location. As the names of the related parameters suggest (InputHeights and OutputHeights), this is often the terrain height. But any other type of information suitable to improve interpolation results is also fine.
Where do we need spatial information soon?
For our OpenDA adapter. Algorithms like the Ensemble Kalman Filter require a spatial noise model. And this model requires the location of certain "parts" of a catchment (e. g. the centres of sub-catchments, hydrological response units or river segments), usually defined within a separate OpenDA XML file. Handling the locations separately from HydPy's project files seems error-prone. Hence, we prefer to implement an HTTP request allowing the adapter to ask the HydPy server for the correct locations.
What do we need to discuss?
- Where is the best place to define location information?
- How much additional support is required?
A very first draw.
We could implement a Location class, allowing us to define x and y and, optionally, z and an epsg code for each Node and Element object. The following (not fully functional) example assumes such a class already exists:
from hydpy import Element, Location, Node
Node("rain_gauge",
location=Location(x=1, y=2, z=3, epsg=25832))
Node("rain_catchment",
location=Location(x=4, y=5, z=6))
Node("outlet_catchment",
location=Location(x=7, y=8))
Element("interpolator",
inlets="rain_gauge",
outlets="rain_catchment")
Element("catchment",
inputs="rain_catchment",
outlets="outlet_catchment",
location=Location(x=9, y=10, z=11))Repeating the same EPSG code multiple time is messy. One way to avoid this would be to provide a default code as an additional option:
pub.options.defaultepsgcode = 25832One could, for example, change the default code at the top of the relevant network file, when required.
In the above example, the "interpolator" element interpolates the precipitation measured at the location of node "rain_gauge" to the location of node "rain_catchment". Node "rain_catchment" passes the interpolated precipitation to element "catchment". Very likely, we want to use the same location information (the area centre of the modelled catchment) for node "rain_gauge" and element "catchment". To not repeatedly define the same location, element "catchment" could query the information from node "rain_gauge" (location="rain_gauge").
Let's assume we want to model the catchment with application model L-Land Version 3. Currently, its control files start like that:
from hydpy.models.lland_v3 import *
parameterstep("1d")
latitude(50.0)
longitude(10.0)
...L-Land Version 3 could query its latitude and its longitude from element "catchment". However, there are two points to consider:
- It would need to transform the given coordinates (in this case UTM) into decimal degrees. Realising this via a library like pyproj does not look to be complicated. However, we should first check if this breaks compatibility with the
arcpylibrary of ArcGIS Pro. Also, we need to make sure it does not complicate our builds or installing HydPy for some other reasons. - Would we convert the control parameters
LongitudeandLatitudeto derived parameters, so that the spatial location never appears within a control file?
The last point connects with the following two problems.
First, if we want OpenDA to modify, for example, the soil moisture of individual hydrological response units, we would ideally also define the location for each unit. Just to write something down:
from hydpy import Element, Location, Node
Element("catchment",
inputs="rain_catchment",
outlets="outlet_catchment",
location=Location(x=9, y=10, z=11),
sublocations=(Location(x=12, y=13), Location(x=14, y=15)))I think this information is too model-specific for a network file and would hamper its readability too much. Hence, if we actually want to support defining separate locations for individual hydrological response units (or river segments...), we should better add a new file type. Maybe, this file type should also support clarifying which variable (e.g. precipitation Nied, soil moisture BoWa, snow water equivalent WAeS) relates to which set of locations, like these two tables:
| LocationType | Variables |
|---|---|
| area centre | Nied |
| hru | BoWa, WAeS |
| Element | LocationType | Coordinates |
|---|---|---|
| catchment | area centre | x=1, y=2, z=3 |
| catchment | hru | (x=4, y=5), (x=6, y=7) |
Second, we also need the locations of selections. Here it seems natural to calculate the average of the locations of all contained elements automatically, weighted with their areas (in case of a subcatchment) or lengths (in case of a river channel). Maybe we could speak of "size" and always assume its value to be given in km² or km:
Element("catchment",
inputs="rain_catchment",
outlets="outlet_catchment",
location=Location(x=9, y=10, z=11),
size=12)| Element | LocationType | Coordinates | Size |
|---|---|---|---|
| catchment | area centre | x=1, y=2, z=3 | 20 |
| catchment | hru | (x=4, y=5), (x=6, y=7) | 5, 15 |
Or, on the Python side, combine the coordinates and the size and also allow defining its unit, for example:
Element("catchment",
inputs="rain_catchment",
outlets="outlet_catchment",
location=Location(x=9, y=10, z=11, size="12km"))
Activity
tyralla commentedon Apr 13, 2021
A few additional thoughts.
The above example regarding the soil moisture (
BoWa) and the snow water equivalent (WAeS) is a little more complicated. HydPy-L-Land requires defining a land-use type for each hydrological response unit. One example isVERS, which stands for sealed areas. There is snow on sealed surfaces but no soil water under them. Hence, we need to take the sealed areas' locations only into account when modifying the snow cover state, not when modifying soil moisture state.A similar problem arises when we want to, for example, calculate the areal mean of soils moisture for a sub-catchment. This is why we allow model developers to define masks. In our example, a mask clarifies for each hydrological response unit if it defines a soil moisture value or not. See the documentation on the method average_values for further information.
Calculating such averages also requires some "size" information. This is not a problem because models like L-Land all implement a parameter (or multiple parameters) defining each hydrological response unit's size. Also, they need to know the total size of the relevant sub-catchment (to convert mm/d to m³/s). One more complete control file example for L-Land:
So, in addition to my initial thoughts, we need to keep in mind that, so far, our model parameter control files are responsible for making spatial information available. And this seems to make sense, as hydrological models represent spatial heterogeneity in many different ways.
Thinking through this perspective, we could add "optional spatial control parameters" (similar to the "optional solver parameters"). For example:
This would allow the model developer to connect different parameters or states to different response units (based on the "mask approach" already available). Also, this would allow dealing with situations where we need to define different locations for different vertical layers (no actual demand so far, but the "vertical structure" of HydPy-W-Land goes into this direction).
Following these thoughts, we could implement a general
HydPy.get_coordinatesmethod, which requires the considered variable types as an argument (and eventually optionally the EPSG code) and returns all relevant coordinates.BGWKlein commentedon Apr 13, 2021
It would be beneficial to add elevation information to the spatial information. This is e.g. requiered for elevation dependent interpolation of temperature. In HBV temperature defined for subbasins is interpolated to the HRUs by using mean elevation of the subbasin and the mean elevation of the HRUs.
tyralla commentedon Apr 13, 2021
Yes, you're right. This is one example, where an available model already handles elevation data on the HRU level (via parameter ZoneZ).
Regarding the subbasin elevation: H-Land comes with three "reference heights". ZRelP is for precipitation, ZRelT for temperature, and ZRelE for evapotranspiration. This allows defining different elevations for the stations where these properties are measured. One further peculiarity of a single model complicating generalisation efforts...
sivogel commentedon Apr 13, 2021
We should keep in mind that a geodata-preprocessing-Toolbox, to read shapefiles/feature classes might be useful in future. We should take this into account for our geospatial information concept.
BGWKlein commentedon Apr 13, 2021
Concerning the different reference heights in HBV: it would be probably better not to use different reference heights for different variables in future and to use only one subbasin reference height provided by the geospatial concept. Otherwise this would over-complicate the geospatial concept.
From my point of view it is better to consider different elevations for the stations in a pre-processing step when interpolating data to the subbasins.
tyralla commentedon Apr 16, 2021
A few notes so as not to forget the key points of our last web meeting:
Element(..., size=...)), we could provide special Element subclasses (e. g.Polygon(..., area=...)),Line(..., length=...)),Point(...))).