Discrepenceis between watershed size reported by RVIC and /watershed API

[corviday]

@sum1lim has found some discrepancies between the number of upstream cells as reported by RVIC and the number reported by the /watershed API.

Lat	Lon	Upstream cells Watershed API	Upstream cells RVIC
-120.8	56.2	1	1
-117.6	55.6	10	10
-121.4	56.1	50	50
-120.3	56	140	2
-121.7	56.4	282	282
-121.6	56.3	360	360
-120.9	56.1	486	486
-120.5	56.2	631	631
-117.7	55.5	1351	2
-117.4	56.1	1948	1941
-122	56	2867	2824
-121.9	56	2868	2752
-121.8	56.1	2900	2893
-121.7	56.1	2925	2918
-121.3	56.2	3300	3293
-121.1	56.2	3342	3903
-120.5	56.1	3921	3335
-120.3	56.1	4558	4550
-120.1	56.1	4708	4700
-118.9	56.1	5104	5096
-118.6	55.9	5248	5140
-118.5	55.9	5256	5247
-118.3	55.9	5264	5256
-118.2	55.9	5266	5258
-118	55.9	5363	5266
-117.8	56	5387	5262
-117.5	56.1	5404	5393

Possibilities:

• domain file has been updated and we're using an out of date version for the watershed API
• issues with numerical parsing, possibly of floating point numbers, they aren't running exactly the same point
• bug in watershed API
• weirdness around grid cells that are partially within a watershed - maybe RVIC ignores this, or doesn't give them full credit? The watershed API counts any square that's in the watershed as a full square.

[corviday]

I have checked to make sure the file that @sum1lim was using for flow direction and the file used by the watershed API for flow direction have the same values. They do.

[corviday]

I'm checking the point at longitude -120.3° latitude 56°. The Watershed API thinks this is a watershed of 140 grid squares; RVIC thinks it is 2 squares.

Here's what the /watershed API thinks the extent of the watershed is. It's generally a reasonably looking watershed: the Kiskatinaw river and tributaries, up to a little bit before it joins the Peace River.

But look at the location of the pin indicating the point that was asked about. Note that this map is so zoomed in the scale is in meters. The selected point is exactly on the boundary.

I think the issue here is that the /watershed API and RVIC are using different rounding algorithms for which square a point is assigned to when it is exactly on the boundary.

Here's the latitude values for the Peace River Flow Direction netCDF:

ncdump -v lat /storage/data/projects/comp_support/climate_explorer_data_prep/hydro/peace_watershed/time-invariant/flow-direction_peace.nc
netcdf flow-direction_peace {
dimensions:
	lon = 192 ;
	lat = 83 ;
variables:
	double lon(lon) ;
	double lat(lat) ;
		lat:standard_name = "latitude" ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
		lat:axis = "Y" ;
	double flow_direction(lat, lon) ;

// global attributes:
    ...

data:

 lat = 53.03125, 53.09375, 53.15625, 53.21875, 53.28125, 53.34375, 53.40625, 
    53.46875, 53.53125, 53.59375, 53.65625, 53.71875, 53.78125, 53.84375, 
    53.90625, 53.96875, 54.03125, 54.09375, 54.15625, 54.21875, 54.28125, 
    54.34375, 54.40625, 54.46875, 54.53125, 54.59375, 54.65625, 54.71875, 
    54.78125, 54.84375, 54.90625, 54.96875, 55.03125, 55.09375, 55.15625, 
    55.21875, 55.28125, 55.34375, 55.40625, 55.46875, 55.53125, 55.59375, 
    55.65625, 55.71875, 55.78125, 55.84375, 55.90625, 55.96875, 56.03125, 
    56.09375, 56.15625, 56.21875, 56.28125, 56.34375, 56.40625, 56.46875, 
    56.53125, 56.59375, 56.65625, 56.71875, 56.78125, 56.84375, 56.90625, 
    56.96875, 57.03125, 57.09375, 57.15625, 57.21875, 57.28125, 57.34375, 
    57.40625, 57.46875, 57.53125, 57.59375, 57.65625, 57.71875, 57.78125, 
    57.84375, 57.90625, 57.96875, 58.03125, 58.09375, 58.15625 ;
}

56 degrees latitude is exactly halfway between 55.96875 and 56.03125. /watershed is assigning this point to the grid cell containing 56.03125, which contains the Kiskatinaw river (140 grid cell watershed). My suspicion is that RVIC is assigning it to the grid cell containing 55.96875, which does not contain a river (2 grid cell watershed).

[sum1lim]

Hmm, I see the point. But there were still some points with a small discrepancy such as (-117.28125, 56.21875) where RVIC thinks it has 7469 upstream grid cells and Watershed API thinks it is 7484. I don't know if it is caused from the same reason.

[corviday]

Hmmmm. It might relate to cells at the edge of the watershed, I think there's some way in a domain file to indicate that a cell is only partly inside the watershed. RVIC might use those fractions to calculate, the area, or it may skip the cells only partly in the watershed. /watershed doesn't do anything like that - either a cell is inside the watershed or it isn't. I'll see if I can turn up anything.

[jameshiebert]

I'd propose that a simple resolution to this issue would consist of not allowing requests on cell boundaries. E.g. if the incoming request coordinate falls on a cell boundary, respond with an HTTP 400 Bad Request error (and an explanatory message).

[corviday]

My first idea about the small differences between RVIC and the watershed API that you mentioned for the whole-watershed value didn't pan out. It would be nice to have more examples. @sum1lim, if you run across any other points that aren't on grid bondaires but don't have the same watershed size, please let me know.

[sum1lim]

@corviday Thanks for investigating on the issue! Following are the locations that show the discrepancies I found until now:

Basin Code | Station ID | Station Name | Lat | Lon | Watershed API | RVIC
BRNFS | 07FC001 | Beatton River near Fort St. John | -120.71875 | 56.28125 | 623 | 622
FRAAR | 07EA005 | Finlay River above Akie River | -125.28125 | 57.15625 | 669 | 666
LSRNG | 07GH002 | Little Smoky River near Guy, AB | -117.15625 | 55.46875 | 433 | 431
ORAOR | 07EC002 | Omineca River above Osilinka River | -124.59375 | 55.90625 | 222 | 219
PERNT | 07FD002 | Peace River near Taylor | -120.65625 | 56.15625 | 3910 | 3836
PRABD | BCGMS | Peace River at Bennett Dam | -122.21875 | 56.03125 | 2836 | 2829
PRAOO | 07EE010 | Pack River at Outlet of McLeod Lake | -123.03125 | 54.96875 | 149 | 148
PRAPR | 07HA001 | Peace River at Peace River, AB | -117.28125 | 56.21875 | 7484 | 7469
SMRAW | 07GJ001 | Smoky River at Watino, AB | -117.59375 | 55.78125 | 1907 | 1900

Lat | Lon | Watershed API | RVIC
-120.3 | 56 | 140 | 2
-117.7 | 55.5 | 1351 | 2
-117.4 | 56.1 | 1948 | 1941
-122 | 56 | 2867 | 2824
-121.9 | 56 | 2868 | 2752
-121.8 | 56.1 | 2900 | 2893
-121.7 | 56.1 | 2925 | 2918
-121.3 | 56.2 | 3300 | 3293
-121.1 | 56.2 | 3342 | 3903
-120.5 | 56.1 | 3921 | 3335
-120.3 | 56.1 | 4558 | 4550
-120.1 | 56.1 | 4708 | 4700
-118.9 | 56.1 | 5104 | 5096
-118.6 | 55.9 | 5248 | 5140
-118.5 | 55.9 | 5256 | 5247
-118.3 | 55.9 | 5264 | 5256
-118.2 | 55.9 | 5266 | 5258
-118 | 55.9 | 5363 | 5266
-117.8 | 56 | 5387 | 5262
-117.5 | 56.1 | 5404 | 5393

[sum1lim]

@corviday , it looks like the small discrepancies only happen on every point with number of upstream grid cells larger than 1000. Points with a number less than 1000 did not show the discrepancies. I cannot confirm this is true, but my observations from the report follow the rule.

Edit: Oh never mind, there are discrepancies between watershed sizes 623 vs 622

Watershed API	RVIC
71	71
67	67
623	622
669	666
84	84
83	83
155	155
361	361
163	163
128	128
104	104
433	431
118	118
90	90
51	51
204	204
179	179
271	271
85	85
222	219
69	69
119	119
3910	3836
2836	2829
444	444
174	174
149	148
7484	7469
47	47
1907	1900
103	103
417	417
41	41

Watershed API	RVIC
1	1
10	10
50	50
282	282
360	360
486	486
631	631
1948	1941
2867	2824
2868	2752
2900	2893
2925	2918
3300	3293
3342	3903
3921	3335
4558	4550
4708	4700
5104	5096
5248	5140
5256	5247
5264	5256
5266	5258
5363	5266
5387	5262
5404	5393