Repository: ec-jrc/lisflood-utilities (language: python)

Issue:
Add utility: water regions map
Add utility to create and verify the water regions map
Add utility: water regions map
Add utility to create and verify the water regions map


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index c123168..1aa4403 100644
--- a/README.md
+++ b/README.md
@@ -40,6 +40,8 @@ Here's a list of utilities you can find in lisflood-utilities package.
 * __compare__ is a package containing a set of simple Python classes that helps to compare 
 netCDF, PCRaster and TSS files.
 
+* __waterregions__ is a package containing two scripts that allow to create and verify a water regions map, respectively.
+
 The package contains convenient classes for reading/writing:
 
 * PCRasterMap
@@ -349,6 +351,109 @@ The Gfit2 tool performs the following steps:
 
 Hosking, J.R.M., 1990. L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. J. R. Stat. Soc. Ser. B Methodol. 52, 105124.
 
+
+## waterregions
+
+The modelling of water abstraction for domestic, industrial, energetic, agricoltural and livestock use   can require a map of the water regions. The concept of water regions and information for their definition are explained [here](htpst://ec-jrc.github.io/lisflood-model/2_18_stdLISFLOOD_water-use/). 
+Since groundwater and surface water resources demand and abstraction are spatially distributed inside each water region, each model set-up must include all the pixels of the water region. This requirement is crucial for the succes of the calibration of the model. This utility allows the user to meet this requirement.
+More specifically, this utility can be used to:
+1. create a water region map which is consistent with a set of calibration points: this purpose is achieved by using the script define_waterregions.
+2. verify the consistency between an existing water region map and an exixting map of calibration catchments: this purpose is achieved by using the script verify_waterregions
+It is here reminded that when calibrating a catchment which is a subset of a larger computational domain, and the option wateruse is switched on, then the option groudwatersmooth must be switched off. The explanation of this requirement is provided in the chapter [Water use](https://ec-jrc.github.io/lisflood-model/2_18_stdLISFLOOD_water-use/) of the LISFLOOD documentation. 
+
+#### Requirements
+python3, pcraster 4.3. The protocol was tested on Linux.
+
+### define_waterregions
+This utility allows to create a  water region map which is consistent with a set of calibration points. The protocol was created by Ad De Roo (Unit D2, Joint Research Centre).
+
+#### Input 
+- List of the coordinates of the calibration points. This list must be provided in a .txt file with three columns: LONGITUDE(or x), LATITUDE(or y), point ID.
+- LDD map in pcraster format. *PCRASTER_VALUESCALE=VS_LDD*. 
+- Countries map in pcraster format. *PCRASTER_VALUESCALE=VS_NOMINAL*. This map shows the political boundaries of the Countries, each Coutry is identified by using a unique ID. This map is used to ensure that the water regions are not split accross different Countries.
+- Map of the initial definition of the water regions in pcraster format. *PCRASTER_VALUESCALE=VS_NOMINAL*. This map is used to attribute a water region to areas not included in the calibration catchments. In order to create this map, the user can follow the guidelines provided [here](https://ec-jrc.github.io/lisflood-model/2_18_stdLISFLOOD_water-use/).
+
+##### Input data provided by this utility:
+This utility provides three maps of [Countries IDs](https://github.com/ec-jrc/lisflood-utilities/tree/master/tests/data): 1arcmin map of Europe (EFAS computational domain), 0.1 degree and 3arcmin maps of the Globe . ACKNOWLEDGEMENTS: both the rasters were retrieved by upsampling the original of the World Borders Datase provided by  http://thematicmapping.org/ (the dataset is available under a Creative Commons Attribution-Share Alike License).
+
+#### Output
+Map of the water regions which is consistent with the calibration catchments. In other words, each water region is entirely included in one calibration catchment.  The test to check the consistency between the newly created water regions map and the calibration catchments is implemented internally by the code and the outcome of the test is printed on the screen. 
+In the output map, each water region is identified by a unique ID.
+
+#### Usage
+The following command line allows to produce a water region map which is consistent with the calibration points:
+
+*python define_waterregions.py -p calib_points_test.txt -l ldd_test.map -C countries_id_test.map -w waterregions_initial_test.map -o my_new_waterregions.map* 
+
+The code internally verifies that the each one of the newly created water regions is entirely included  within one calibration catchments. If this condition is satisfied, the follwing message in printed out: *“OK! Each water region is completely included inside one calibration catchment”*. If the condition is not satisfied, the error message is *“ERROR: The  water regions WR are included in more than one calibration catchment”*. Moreover, the code provides the list of the water regions WR and the calibration catchments that do not meet the requirment. This error highlight a problem in the input data: the user is recommended to check (and correct) the list of calibration points and the input maps.
+
+The input and output arguments are listed below. All the inputs are required. 
+
+
+```text
+usage: define_waterregions.py [-h] -p CALIB_POINTS -l LDD -C COUNTRIES_ID -w
+                              WATERREGIONS_INITIAL -o OUTPUT_WR
+
+Define Water Regions consistent with calibration points: {}
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -p CALIB_POINTS, --calib_points CALIB_POINTS
+                        list of calibration points: lon or x, lat or y, point id
+  -l LDD, --ldd LDD     LDD map pcraster format
+  -C COUNTRIES_ID, --countries_id COUNTRIES_ID
+                        map of Countries ID pcraster format
+  -w WATERREGIONS_INITIAL, --waterregions_initial WATERREGIONS_INITIAL
+                        initial map of water regions pcraster format
+  -o OUTPUT_WR, --output_wr OUTPUT_WR
+                        output map of water regions pcraster format
+```
+
+NOTE: 
+The utility **nc2pcr** can be used to convert a netcdf map into pcraster format.
+
+
+### verify_waterregions
+
+This function allows to verify the consistency between a  water region map and a map of calibration catchments. This function must be used when the water region map and the map of calibration catchments have been defined in an independent manner (i.e. not using the utility **define_waterregions**). The function verify_waterregions verifies that each water region map is entirely included in one calibration catchment. If this condition is not satisfied, an error message is printed on the screen. 
+
+#### Input
+- Map of calibration catchments in netcdf format.
+- Water regions map in netcdf format.
+
+#### Output
+The output is a message on the screen. There are two options:
+- 'OK! Each water region is completely included inside one calibration catchment.'
+- 'ERROR: The  water regions WR are included in more than one calibration catchment’: this message is followed by the list of the water regions and of the catchment that raised the isuue.
+In case of error message, the user can implement the function **define_waterregions**.
+
+#### Usage
+The following command line allows to produce a water region map which is consistent with the calibration points:
+
+*python verify_waterregions.py -cc calib_catchments_test.nc -wr waterregions_test.nc*
+
+The input and output arguments are listed below. All the inputs are required. 
+
+```text
+usage: verify_waterregions.py [-h] -cc CALIB_CATCHMENTS -wr WATERREGIONS
+
+Verify that the Water Regions map is consistent with the map of the
+calibration catchments
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -cc CALIB_CATCHMENTS, --calib_catchments CALIB_CATCHMENTS
+                        map of calibration catchments, netcdf format
+  -wr WATERREGIONS, --waterregions WATERREGIONS
+                        map of water regions, netcdf format
+```
+
+NOTE:
+The utility **pcr2nc** can be used to convert a map in pcraster format into netcdf format.
+
+
+
+
 ## Using lisfloodutilities programmatically TODO
 
 You can use lisflood utilities in your python programs.
diff --git a/bin/define_waterregions b/bin/define_waterregions
new file mode 100755
index 0000000..b5ebc9a
--- /dev/null
+++ b/bin/define_waterregions
@@ -0,0 +1,14 @@
+#!python
+
+import os
+import sys
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+src_dir = os.path.normpath(os.path.join(current_dir, '../src/'))
+if os.path.exists(src_dir):
+    sys.path.append(src_dir)
+
+from lisfloodutilities.waterregions.define_waterregions import main_script
+
+if __name__ == '__main__':
+    main_script()
diff --git a/bin/verify_waterregions b/bin/verify_waterregions
new file mode 100755
index 0000000..310fc23
--- /dev/null
+++ b/bin/verify_waterregions
@@ -0,0 +1,14 @@
+#!python
+
+import os
+import sys
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+src_dir = os.path.normpath(os.path.join(current_dir, '../src/'))
+if os.path.exists(src_dir):
+    sys.path.append(src_dir)
+
+from lisfloodutilities.waterregions.very_waterregions import main_script
+
+if __name__ == '__main__':
+    main_script()
diff --git a/src/lisfloodutilities/waterregions/__init__.py b/src/lisfloodutilities/waterregions/__init__.py
new file mode 100644
index 0000000..90c039e
--- /dev/null
+++ b/src/lisfloodutilities/waterregions/__init__.py
@@ -0,0 +1,17 @@
+"""
+
+Copyright 2019-2020 European Union
+
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+
diff --git a/src/lisfloodutilities/waterregions/define_waterregions.py b/src/lisfloodutilities/waterregions/define_waterregions.py
new file mode 100755
index 0000000..6f56913
--- /dev/null
+++ b/src/lisfloodutilities/waterregions/define_waterregions.py
@@ -0,0 +1,138 @@
+"""
+Copyright 2019-2020 European Union
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+
+import argparse
+import sys
+import numpy as np
+import pcraster as pcr
+import os
+import time
+import logging
+
+logging.basicConfig(format='[%(asctime)s] - %(message)s', datefmt='%H:%M:%S', level=logging.INFO)
+logger = logging.getLogger()
+
+
+def main(cliargs):
+    parser = ParserHelpOnError(description='Define Water Regions consistent with calibration points: {}')
+    parser.add_args()
+    args = parser.parse_args(cliargs)
+
+    calib_points = args.calib_points
+    countries_id = args.countries_id
+    ldd = args.ldd
+    waterregions_initial = args.waterregions_initial
+    output_wr = args.output_wr
+
+    [waterregion, subcat1] = define_waterregions(calib_points, countries_id, ldd, waterregions_initial,output_wr)
+    logger.info('\nUsing %s and %s to define the water regions\n ', calib_points, ldd)
+    
+    # check the consistency between the water regions and the calibration catchments
+    cal_catch = pcr.pcr_as_numpy(subcat1)
+    cal_catch[cal_catch < 1] = -9999
+    cal_catch[np.isnan(cal_catch) == 1] = -9999
+    
+    wr = pcr.pcr_as_numpy(waterregion)
+    wr_id = np.unique(wr)
+       
+    id_error_wr = []
+    cal_catch_error_wr = []
+    output_message = []
+  
+    for a in wr_id:
+           extract_wr = np.where(wr == a,cal_catch,-9999)
+           num_cal_catch = np.unique(extract_wr)
+           num_cal_catch_check = np.extract(num_cal_catch != -9999,num_cal_catch)
+           if len(num_cal_catch_check) > 1:
+             a = int(a)
+             id_error_wr.append(a)
+             num_cal_catch_write = num_cal_catch[num_cal_catch != -9999].astype(int)
+             cal_catch_error_wr.append(num_cal_catch_write)
+             output_message = 'ERROR: The  water regions WR are included in more than one calibration catchment \n WR=' + str(id_error_wr) + '\n Calibration catchments= ' + str(cal_catch_error_wr)
+    if id_error_wr == []:
+             output_message = 'OK! Each water region is completely included inside one calibration catchment.'
+    
+    print(output_message)
+    os.unlink('./points.map')
+
+
+class ParserHelpOnError(argparse.ArgumentParser):
+    def error(self, message):
+        sys.stderr.write('Error: %s\n' % message)
+        self.print_help()
+        sys.exit(1)
+        
+    def add_args(self):
+
+        self.add_argument("-p", "--calib_points",
+                          help='list of calibration points: lon or x, lat or y, point id',
+                          required=True)
+        self.add_argument("-l", "--ldd",
+                          help='LDD map pcraster format',
+                          required=True)
+        self.add_argument("-C", "--countries_id",
+                          help='map of Countries ID pcraster format',
+                          required=True)
+        self.add_argument("-w", "--waterregions_initial",
+                          help='initial map of water regions pcraster format',
+                          required=True)
+        self.add_argument("-o", "--output_wr",
+                          help='output map of water regions pcraster format',
+                          required=True)
+
+def define_waterregions(calib_points=None, countries_id=None, ldd=None, waterregions_initial=None, output_wr=None):
+
+    
+    #1. The calibration points are converted into a map
+    command_string = 'col2map -N ' + calib_points + ' points.map --large --clone ' + countries_id
+    os.system(command_string)
+
+    #2. Cacthment map1 derived from calibration points
+    ldd1 = pcr.readmap(ldd)
+    points = pcr.readmap('points.map')
+    subcat1 = pcr.subcatchment(ldd1,points)
+   
+    #3. Making map with all valid ldd cells
+    land = pcr.scalar(subcat1)
+    land = pcr.nominal(land >= 0.99999)
+  
+    #4. Taking only the catchments from calibration points
+    subcat2 = pcr.ifthen(subcat1 != 0,subcat1) 
+
+    #5. Splitting riverbasins if part of different countries, so if cross-border catchments
+    scalar_a = pcr.scalar(subcat2)
+    scalar_b = pcr.scalar(countries_id)
+    subcat3 = pcr.nominal(scalar_a*scalar_b)
+
+    #6. Covering the missing areas with the old waterregions
+    subcat4 = pcr.cover(subcat3,waterregions_initial)
+
+    #7. Make sure that all land cells are filled
+    subcat5 = pcr.cover(subcat4,land)
+    
+    #8. Each water region is given a unique ID number
+    waterregion0 = pcr.nominal(subcat5)
+    waterregion = pcr.clump(waterregion0)
+
+
+    #9. Save the water region map
+    pcr.report(waterregion,output_wr)
+    
+    return waterregion, subcat1
+                          
+def main_script():
+    sys.exit(main(sys.argv[1:]))
+
+
+if __name__ == '__main__':
+    main_script()                         
+
diff --git a/src/lisfloodutilities/waterregions/verify_waterregions.py b/src/lisfloodutilities/waterregions/verify_waterregions.py
new file mode 100755
index 0000000..374e544
--- /dev/null
+++ b/src/lisfloodutilities/waterregions/verify_waterregions.py
@@ -0,0 +1,97 @@
+"""
+Copyright 2019-2020 European Union
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+
+import argparse
+import sys
+import numpy as np
+import os
+import time
+import logging
+from netCDF4 import Dataset
+
+logging.basicConfig(format='[%(asctime)s] - %(message)s', datefmt='%H:%M:%S', level=logging.INFO)
+logger = logging.getLogger()
+
+
+def main(cliargs):
+    parser = ParserHelpOnError(description='Verify that the Water Regions map is consistent with the map of the calibration catchments')
+    parser.add_args()
+    args = parser.parse_args(cliargs)
+    calib_catchments = args.calib_catchments 
+    waterregions = args.waterregions 
+    output_message = verify_waterregions(calib_catchments, waterregions)
+    logger.info('\n Verify the consistency between %s and %s ', calib_catchments, waterregions)
+    logger.info(output_message)
+
+
+class ParserHelpOnError(argparse.ArgumentParser):
+    def error(self, message):
+        sys.stderr.write('Error: %s\n' % message)
+        self.print_help()
+        sys.exit(1)
+        
+    def add_args(self):
+        self.add_argument("-cc", "--calib_catchments",
+                          help='map of calibration catchments, netcdf format',
+                          required=True)
+        self.add_argument("-wr", "--waterregions",
+                          help='map of water regions, netcdf format',
+                          required=True)
+
+
+def verify_waterregions(calib_catchments=None, waterregions=None):         
+       waterregions_map = Dataset(waterregions,'r','format=NETCDF4_classic')
+       
+       for v in waterregions_map.variables:
+           if len(waterregions_map.variables[v].dimensions) == 2:
+              wr = waterregions_map.variables[v][:]             
+      
+       calibration_catchments_map = Dataset(calib_catchments,'r','format=NETCDF4_classic')
+       for v in calibration_catchments_map.variables:
+           if len(calibration_catchments_map.variables[v].dimensions) == 2:
+              cal_catch = calibration_catchments_map.variables[v][:]
+       
+       
+       cal_catch[cal_catch<1] = -9999
+       cal_catch[np.isnan(cal_catch) == 1] = -9999
+       wr_id = np.unique(wr)
+       
+       id_error_wr = []
+       cal_catch_error_wr = []
+       output_message = []
+  
+       for a in wr_id:
+           extract_wr = np.where(wr == a, cal_catch, -9999)
+           num_cal_catch = np.unique(extract_wr)
+           num_cal_catch_check = np.extract(num_cal_catch != -9999, num_cal_catch)
+           if len(num_cal_catch_check) > 1:
+             id_error_wr.append(a)
+             num_cal_catch_write = num_cal_catch[num_cal_catch != -9999].astype(int)
+             cal_catch_error_wr.append(num_cal_catch_write)
+             output_message = 'ERROR: The  water regions WR are included in more than one calibration catchment \n WR=' + str(id_error_wr) + '\n Calibration catchments= ' + str(cal_catch_error_wr)
+       if id_error_wr == []:
+           output_message = 'PASSED: Each calibration catchment contains only a finite number of water regions.'
+    
+
+       waterregions_map.close()
+       calibration_catchments_map.close()
+       
+ 
+       return output_message
+                          
+def main_script():
+    sys.exit(main(sys.argv[1:]))
+
+
+if __name__ == '__main__':
+    main_script()                         
+


Test patch:
diff --git a/tests/data/waterregions/Europe_countries_id_1arcmin.nc b/tests/data/waterregions/Europe_countries_id_1arcmin.nc
new file mode 100755
index 0000000..73e0bf4
Binary files /dev/null and b/tests/data/waterregions/Europe_countries_id_1arcmin.nc differ
diff --git a/tests/data/waterregions/Global_countries_id_01deg.nc b/tests/data/waterregions/Global_countries_id_01deg.nc
new file mode 100755
index 0000000..a712d5c
Binary files /dev/null and b/tests/data/waterregions/Global_countries_id_01deg.nc differ
diff --git a/tests/data/waterregions/Global_countries_id_3arcmin.nc b/tests/data/waterregions/Global_countries_id_3arcmin.nc
new file mode 100755
index 0000000..d58f925
Binary files /dev/null and b/tests/data/waterregions/Global_countries_id_3arcmin.nc differ
diff --git a/tests/data/waterregions/calib_catchments_test.nc b/tests/data/waterregions/calib_catchments_test.nc
new file mode 100755
index 0000000..cc6157a
Binary files /dev/null and b/tests/data/waterregions/calib_catchments_test.nc differ
diff --git a/tests/data/waterregions/calib_points_test.txt b/tests/data/waterregions/calib_points_test.txt
new file mode 100755
index 0000000..d437fc5
--- /dev/null
+++ b/tests/data/waterregions/calib_points_test.txt
@@ -0,0 +1,245 @@
+-80.65 -5.25 1
+-80.25 -4.35 2
+-79.45 0.45 3
+-78.95 9.15 4
+-77.55 -4.45 5
+-76.65 3.05 6
+-76.55 3.15 7
+-76.55 -6.45 8
+-76.45 3.25 9
+-76.45 3.35 10
+-76.35 3.85 11
+-76.25 -6.55 12
+-76.15 4.35 13
+-75.95 4.85 14
+-75.85 6.45 15
+-75.75 5.85 16
+-75.65 2.25 17
+-75.65 5.25 18
+-75.55 2.45 19
+-75.35 2.95 20
+-75.35 7.35 21
+-75.15 3.35 22
+-74.95 3.75 23
+-74.95 8.05 24
+-74.85 4.25 25
+-74.75 5.15 26
+-74.75 8.15 27
+-74.75 9.25 28
+-74.45 6.35 29
+-74.05 6.85 30
+-73.95 7.35 31
+-73.95 8.95 32
+-73.95 -4.45 33
+-73.85 8.65 34
+-73.65 -4.75 35
+-73.15 -16.35 36
+-73.15 -4.05 37
+-73.05 -3.45 38
+-72.85 -5.15 39
+-72.75 -8.95 40
+-72.65 -7.65 41
+-72.45 -13.25 42
+-72.45 -15.95 43
+-71.75 -7.05 44
+-71.45 -4.65 45
+-70.95 -4.35 46
+-70.75 -8.15 47
+-70.45 -8.25 48
+-70.45 -4.95 49
+-70.25 -4.55 50
+-70.05 -7.45 51
+-69.95 -15.25 52
+-69.95 -6.75 53
+-69.95 -4.25 54
+-69.45 -1.35 55
+-69.35 -15.95 56
+-69.15 -5.55 57
+-69.15 -8.75 58
+-69.15 0.45 59
+-68.85 -9.75 60
+-68.75 -3.45 61
+-68.55 -9.05 62
+-68.35 -6.55 63
+-68.15 -4.75 64
+-67.95 -4.25 65
+-67.95 -3.15 66
+-67.85 -10.05 67
+-67.65 1.15 68
+-67.55 -7.55 69
+-67.55 -3.35 70
+-67.55 -14.55 71
+-67.35 -4.35 72
+-67.35 -9.05 73
+-67.35 -8.65 74
+-67.35 0.45 75
+-67.25 -11.45 76
+-66.95 -7.75 77
+-66.85 1.35 78
+-66.85 -4.95 79
+-66.85 -0.15 80
+-66.75 -11.65 81
+-66.65 -1.85 82
+-66.05 -10.95 83
+-65.55 -9.85 84
+-65.55 -10.55 85
+-65.35 -10.75 86
+-65.35 -9.75 87
+-65.25 -4.35 88
+-65.05 -13.05 89
+-64.95 -6.35 90
+-64.95 -15.35 91
+-64.95 -0.45 92
+-64.75 -7.25 93
+-63.75 -8.75 94
+-63.15 3.45 95
+-63.05 -7.45 96
+-63.05 -4.05 97
+-62.95 -12.85 98
+-62.25 -4.75 99
+-62.05 -8.95 100
+-61.95 -10.85 101
+-61.75 -6.55 102
+-61.45 2.85 103
+-61.35 -5.85 104
+-61.25 -11.75 105
+-61.15 -11.65 106
+-61.15 3.45 107
+-61.05 -13.55 108
+-61.05 1.85 109
+-60.95 2.55 110
+-60.75 -7.65 111
+-60.65 -9.85 112
+-60.65 -7.25 113
+-60.65 -3.35 114
+-60.55 -7.85 115
+-60.55 3.25 116
+-60.55 -7.75 117
+-60.35 6.35 118
+-59.95 -15.05 119
+-59.65 -3.05 120
+-59.45 -10.25 121
+-59.45 -1.95 122
+-59.05 -6.85 123
+-59.05 -6.65 124
+-58.75 5.35 125
+-58.55 -1.55 126
+-58.35 -11.45 127
+-58.25 -7.55 128
+-58.15 -7.35 129
+-58.05 -0.65 130
+-58.05 -11.75 131
+-57.95 -10.65 132
+-57.85 -1.45 133
+-57.75 -10.85 134
+-57.45 -11.55 135
+-57.45 -8.85 136
+-57.45 -5.85 137
+-57.15 -1.15 138
+-57.05 -1.05 139
+-56.85 -0.75 140
+-56.15 -0.95 141
+-56.15 -9.65 142
+-55.95 -4.25 143
+-55.85 -12.65 144
+-55.75 -6.25 145
+-55.65 -11.55 146
+-55.65 -10.95 147
+-55.55 -10.15 148
+-55.55 -1.95 149
+-55.35 -13.55 150
+-55.25 -10.15 151
+-54.95 -1.45 152
+-54.55 -5.85 153
+-54.45 4.95 154
+-53.95 -4.45 155
+-53.05 -13.45 156
+-53.05 -15.25 157
+-52.55 -0.55 158
+-52.35 -14.65 159
+-52.25 -6.95 160
+-52.15 0.95 161
+-52.15 -3.35 162
+-51.85 3.85 163
+-51.75 -14.05 164
+-51.75 -6.75 165
+-51.55 -3.65 166
+-51.55 0.75 167
+-51.05 -14.85 168
+-50.95 -12.25 169
+-50.75 -13.65 170
+-50.65 -11.85 171
+-50.65 -11.65 172
+-50.55 -13.25 173
+-50.35 -13.55 174
+-49.95 -12.05 175
+-49.95 -10.55 176
+-49.55 -12.45 177
+-49.55 -15.25 178
+-49.25 -5.45 179
+-49.25 -8.35 180
+-49.05 -14.75 181
+-49.05 -14.55 182
+-49.05 -5.35 183
+-48.65 -14.95 184
+-48.65 -11.95 185
+-48.55 -12.05 186
+-48.55 -6.35 187
+-48.35 -9.55 188
+-48.25 -11.55 189
+-48.25 -12.75 190
+-48.15 -8.35 191
+-47.95 -9.35 192
+-47.85 -21.05 193
+-47.75 -2.55 194
+-47.65 -9.95 195
+-47.45 -7.35 196
+-47.35 -7.75 197
+-46.95 -3.45 198
+-46.85 -13.75 199
+-46.25 -1.95 200
+-45.65 -8.25 201
+-45.45 -3.65 202
+-45.35 -16.95 203
+-45.25 -19.25 204
+-45.25 -3.75 205
+-45.15 -16.15 206
+-45.05 -16.65 207
+-45.05 -16.35 208
+-44.95 -17.25 209
+-44.85 -7.05 210
+-44.75 -7.35 211
+-44.65 -7.45 212
+-44.65 -4.55 213
+-44.25 -3.65 214
+-44.15 -14.25 215
+-44.05 -11.55 216
+-43.95 -3.45 217
+-43.95 -11.25 218
+-43.95 -14.75 219
+-43.75 -14.25 220
+-43.35 -4.85 221
+-43.25 -11.55 222
+-42.95 -6.75 223
+-42.85 -5.15 224
+-42.55 -5.45 225
+-42.25 -3.95 226
+-42.25 -5.65 227
+-42.25 -3.45 228
+-42.15 -5.75 229
+-41.95 -3.65 230
+-41.75 -5.25 231
+-41.65 -18.95 232
+-40.55 -9.45 233
+-40.35 -16.15 234
+-39.85 -12.55 235
+-39.65 -10.95 236
+-39.45 -12.65 237
+-39.35 -15.55 238
+-39.25 -12.05 239
+-38.85 -6.45 240
+-38.35 -5.15 241
+-37.85 -11.75 242
+-37.35 -6.45 243
+-37.25 -11.05 244
+-34.95 -7.15 245
diff --git a/tests/data/waterregions/countries_id_test.map b/tests/data/waterregions/countries_id_test.map
new file mode 100755
index 0000000..7e51a8f
Binary files /dev/null and b/tests/data/waterregions/countries_id_test.map differ
diff --git a/tests/data/waterregions/ldd_test.map b/tests/data/waterregions/ldd_test.map
new file mode 100755
index 0000000..06b8b0d
Binary files /dev/null and b/tests/data/waterregions/ldd_test.map differ
diff --git a/tests/data/waterregions/waterregions_initial_test.map b/tests/data/waterregions/waterregions_initial_test.map
new file mode 100755
index 0000000..e839f05
Binary files /dev/null and b/tests/data/waterregions/waterregions_initial_test.map differ
diff --git a/tests/data/waterregions/waterregions_reference_output.map b/tests/data/waterregions/waterregions_reference_output.map
new file mode 100755
index 0000000..0bfbded
Binary files /dev/null and b/tests/data/waterregions/waterregions_reference_output.map differ
diff --git a/tests/data/waterregions/waterregions_test.nc b/tests/data/waterregions/waterregions_test.nc
new file mode 100755
index 0000000..4fd9473
Binary files /dev/null and b/tests/data/waterregions/waterregions_test.nc differ
diff --git a/tests/data/waterregions/waterregions_test_NOTc.nc b/tests/data/waterregions/waterregions_test_NOTc.nc
new file mode 100644
index 0000000..ecc3b3a
Binary files /dev/null and b/tests/data/waterregions/waterregions_test_NOTc.nc differ
diff --git a/tests/test_waterregions.py b/tests/test_waterregions.py
new file mode 100755
index 0000000..9591d6e
--- /dev/null
+++ b/tests/test_waterregions.py
@@ -0,0 +1,41 @@
+import os
+
+from lisfloodutilities.readers import PCRasterMap
+
+from lisfloodutilities.waterregions.define_waterregions import define_waterregions
+from lisfloodutilities.waterregions.verify_waterregions import verify_waterregions
+
+import pcraster as pcr
+
+class TestWaterRegions:
+
+    def test_define_waterregions(self):
+        calib_points = 'tests/data/waterregions/calib_points_test.txt'
+        countries_id = 'tests/data/waterregions/countries_id_test.map'
+        ldd = 'tests/data/waterregions/ldd_test.map'
+        waterregion_init = 'tests/data/waterregions/waterregions_initial_test.map'
+        output_wr = 'tests/data/waterregions/waterregions_output_test.map'
+        reference_output = 'tests/data/waterregions/waterregions_reference_output.map'
+
+        define_waterregions(calib_points, countries_id, ldd, waterregion_init, output_wr)
+
+        out = PCRasterMap(output_wr)
+        reference = PCRasterMap(reference_output)
+
+        assert out == reference
+
+        out.close()
+        reference.close()
+        os.unlink(output_wr)
+        
+    def test_verify_waterregions(self):
+        calib_catchments = 'tests/data/waterregions/calib_catchments_test.nc'
+        waterregions = 'tests/data/waterregions/waterregions_test.nc'
+        output_message = verify_waterregions(calib_catchments, waterregions)
+        assert output_message == 'PASSED: Each calibration catchment contains only a finite number of water regions.'
+
+    def test_verify_waterregions_fail(self):
+        calib_catchments = 'tests/data/waterregions/calib_catchments_test.nc'
+        waterregions = 'tests/data/waterregions/waterregions_test_NOTc.nc'
+        output_message = verify_waterregions(calib_catchments, waterregions)
+        assert output_message != 'PASSED: Each calibration catchment contains only a finite number of water regions.'