widths seem off for rivers that start outside the model domain and for reaches with multiple upstream tributaries

[janetrbarclay]

I am using sfrmaker version 0.7.0 and am working a model with multiple rivers that begin upstream of the model domain. After creating the lines instances from NHDPlus, when I run lns.to_sfr(), it appears the asum is reset to 0 at the model boundary, regardless of any out-of-boundary length. In addition, it seems that only 1 upstream reach is included in calculating asum1, even if there should be multiples.

For example, for line_id = 7701178 (the CT River in southern CT), asum2 in the lines instance is 18782002 meters, but in the sfrdata instance it is reduced to just a few hundred meters (with significant implications for the calculated widths).

lns = sfrmaker.Lines.from_nhdplus_v2(NHDPlus_paths=NHDPlus_paths, filter =os.path.join(geo_ws,"clip_box.shp")) lns.df.loc[lns.df.id==7701178] Out[19]: id toid asum1 asum2 7701178 7701178 [7701204] 0 18782002.0

sfrdata = lns.to_sfr(model=gwf, model_length_units='meters', active_area=os.path.join(geo_ws,"domain_outline.shp")) sfrdata.reach_data.loc[sfrdata.reach_data.line_id==7701178,['line_id','ireach','rchlen','width','asum']] Out[17]: line_id ireach rchlen width asum 1727 7701178 1 18.881460 1.0 9.440730 1728 7701178 2 354.118652 1.0 195.940781 1729 7701178 3 38.035713 1.0 392.017944

I can calculate the widths for the lines instance explicitly using the example in lines.to_sfr(), but in the existing code it seems that only 1 upstream segment is included in the routing dictionary, even if there are multiple. Again using line_id = 7701178 as an example:

routing_r = {v:k for k,v in lns.routing.items() if v!=0} routing_r[7701178]

Picking routing connections at divergences... finished in 0.21s

Out[21]: 7701126

But there should be two upstream reaches (and the one that's missing is the more important one for calculating asum1 to the reach): id toid asum1 asum2
7701124 7701124 [7701178] 0 18769223.0 7701126 7701126 [7701178] 0 8163.0

In case it's helpful, here's the alternative that I'm using to include all the upstream reaches:

routing_r = {v:[] for k,v in lns.routing.items() if v!=0}
for k1,v1 in lns.routing.items():
   if v1!=0:
       routing_r[v1].append(k1)

asums = dict(zip(lns.df.id,lns.df.asum2 * sfrmaker.units.convert_length_units(lns.attr_length_units,'meters')))
lns.df['asum1'] = [np.sum([asums.get(x, 0) for x in routing_r.get(id, 0)]) if routing_r.get(id, 0)!=0 else 0 for id in lns.df.id.values]

Using this I get both upstream reaches in the routing dictionary: routing_r[7701178] Out[34]: [7701124, 7701126]

and asum1 includes both reaches: lns.df.loc[lns.df.id==7701178]

Out[35]: id toid asum1 asum2 width1 width2 \ 7701178 7701178 [7701204] 18777386.0 18782002.0 166.245005 166.265568

[aleaf]

Hi @janetrbarclay, sorry for the slow reply. Thanks for bringing this to my attention! I think this issue should be fixed with 2c8978b1. Can you confirm? Thanks

[janetrbarclay]

No worries, @aleaf. I have a full schedule today and tomorrow, but will check it on Friday. Thanks!

[janetrbarclay]

HI @aleaf. Yes, it looks like it is fixed. Now when create a lines instance and check that reach I get this (as expected):

lns = sfrmaker.Lines.from_nhdplus_v2(NHDPlus_paths=NHDPlus_paths, filter =os.path.join(geo_ws,"clip_box.shp")) lns.df.loc[lns.df.id==7701178]

Out[24]: id toid asum1 asum2 width1 width2 elevup elevdn name geometry 7701178 7701178 [7701204] 18777385.0 18782002.0 0 0 0.0 0.0 Connecticut River LINESTRING Z (-72.62282279771489 41.6012024020...

and then after sfrdata = lns.to_sfr(model=gwf, model_length_units='meters', active_area=os.path.join(geo_ws,"domain_outline.shp")) sfrdata.reach_data.loc[sfrdata.reach_data.line_id==7701178,['line_id','ireach','rchlen','width','asum']]

Out[18]: line_id ireach rchlen width asum 1740 7701178 1 117.990074 166.245270 18777444.0 1741 7701178 2 293.045746 166.246185 18777650.0

which is also as expected. Thanks so much for your help!