search

MiKatt / openSTARS

open source implementation of the STARS ArcGIS toolbox
https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0239237
Other
46 stars 13 forks source link

Multiple stream identifiers fix for calcCatchmArea_assignNetID #27

Closed sigmafelix closed 2 years ago

sigmafelix commented 2 years ago

Hello, thank you for maintaining the package. I found that there is an error when the calc_edges function tries to assign network IDs using calcCatchmArea_assignNetID.

calcCatchmArea_assignNetID = 
function (dt, id, net_ID) 
{
    **if (dt[stream == id, prev_str01, ] == 0)** {
        dt[stream == id, `:=`(total_area, area)]
        dt[stream == id, `:=`(netID, net_ID)]
    }
    else {
        a1 <- calcCatchmArea_assignNetID(dt, dt[stream == id, 
            prev_str01], net_ID)
        a2 <- calcCatchmArea_assignNetID(dt, dt[stream == id, 
            prev_str02], net_ID)
        dt[stream == id, `:=`(total_area, a1 + a2 + dt[stream == 
            id, area])]
        dt[stream == id, `:=`(netID, net_ID)]
    }
    return(dt[stream == id, total_area])
}

An error the condition has length > 1 and only the first element will be used is printed at the line highlighted with the double asterisks in the code block above when the dt.streams object happens to have multiple rows with the same stream value. I guess the highlighted line should be changed to have if (sum(dt[stream == id, prev_str01, ]) == 0) to cover the multiple-stream values cases.

Thank you very much.

MiKatt commented 2 years ago

Thanks for bringing this to my attention. Actually, stream should be unique, so the condition should always have length 1. Could you provide me a reproducible example? If the condition has length > 1, there is something wrong in the functions elsewhere and I would like to fix it.

sigmafelix commented 2 years ago

Thank you for your prompt response! Here is my reprex code:

download.file('https://github.com/sigmafelix/Gmisc/raw/master/AZ_highcomp.tif', 'AZ_highcomp.tif')
download.file('https://github.com/sigmafelix/Gmisc/raw/master/sites_az_reprojected.geojson', 'sites_az_reprojected.geojson')

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.2
#> WARNING: different compile-time and runtime versions for GEOS found:
#> Linked against: 3.8.0-CAPI-1.13.1  compiled against: 3.8.1-CAPI-1.13.3
#> It is probably a good idea to reinstall sf, and maybe rgeos and rgdal too
library(raster)
#> Loading required package: sp
library(openSTARS)
#> Loading required package: data.table
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:raster':
#> 
#>     shift
#> Loading required package: rgrass7
#> Loading required package: XML
#> GRASS GIS interface loaded with GRASS version: GRASS 7.8.2 (2019)
#> and location: file890e166c4d1eb
#> Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.8.0-CAPI-1.13.1
#> and GEOS at installation 3.8.1-CAPI-1.13.3differ
library(rgrass7)

setup_grass_environment(dem = 'AZ_highcomp.tif', gisBase = '/usr/lib/grass78/', override = TRUE)
import_data(dem = 'AZ_highcomp.tif', sites = 'sites_az_reprojected.geojson')
#> Loading DEM into GRASS as 'dem' ...
#> Loading sites into GRASS as 'sites_o' ...
#> No streams available, skipping.
derive_streams(accum_threshold = 100)
#> Conditioning DEM ...
#> Load elevation map
#>    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
#> Remove one cell extremas
#>    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
#> Set edge points
#>    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
#> A* Search...
#>    0%   2%   4%   6%   8%  10%  12%  14%  16%  18%  20%  22%  24%  26%  28%  30%  32%  34%  36%  38%  40%  42%  44%  46%  48%  50%  52%  54%  56%  58%  60%  62%  64%  66%  68%  70%  72%  74%  76%  78%  80%  82%  84%  86%  88%  90%  92%  94%  96%  98% 100%
#> Processing 14320 sinks
#>    1%   3%   5%   7%   9%  11%  13%  15%  17%  19%  21%  23%  25%  27%  29%  31%  33%  35%  37%  39%  41%  43%  45%  47%  49%  51%  53%  55%  57%  59%  61%  63%  65%  67%  69%  71%  73%  75%  77%  79%  81%  83%  85%  87%  89%  91%  93%  95%  97%  99% 100%
#> Write conditioned DEM
#>    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
#> Deriving streams from DEM ...
#> Calculating stream topology ...
#> WARNING: Stream network may be too dense
#> WARNING: Stream network may be too dense
#> Derived streams saved as 'streams_v'.
correct_compl_confluences()
#> Fixing 2 complex confluences with 4 upstream segments each
#> Original stream topology file moved to 'streams_v_o4'.
#> Breaking lines and moving vertices ...
#> WARNING: Values in column <length_new> will be overwritten
#> Updating topology ...
#> Fixing 115 complex confluences with 3 upstream segments each
#> Original stream topology file moved to 'streams_v_o3'.
#> Breaking lines and moving vertices ...
#> WARNING: Values in column <length_new> will be overwritten
#> Updating topology ...
#> Complex confluences were removed. Please check changed features in 'streams_v'.
calc_edges()
#> Calculating reach contributing area (RCA) ...
#> Calculating upstream catchment areas ...
#> Warning in if (dt[stream == id, prev_str01, ] == 0) {: the condition has length
#> > 1 and only the first element will be used
#> Warning in if (dt[stream == id, prev_str01, ] == 0) {: the condition has length
#> > 1 and only the first element will be used

#> Warning in if (dt[stream == id, prev_str01, ] == 0) {: the condition has length
#> > 1 and only the first element will be used

#> Warning in if (dt[stream == id, prev_str01, ] == 0) {: the condition has length
#> > 1 and only the first element will be used
#> Error in .prepareFastSubset(isub = isub, x = x, enclos = parent.frame(), : RHS of == is length 2 which is not 1 or nrow (17310). For robustness, no recycling is allowed (other than of length 1 RHS). Consider %in% instead.
#calc_sites()
MiKatt commented 2 years ago

Would you please update all R packages and check if the error still persists?

sigmafelix commented 2 years ago

I tried to update all packages I have and still had the same issue. here is my sessionInfo() result:

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] openSTARS_1.2.2   rgrass7_0.2-6     XML_3.99-0.8      data.table_1.14.2
[5] raster_3.5-2      sp_1.4-5          sf_1.0-3         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.7         pillar_1.6.4       compiler_4.1.1     prettyunits_1.1.1 
 [5] progress_1.2.2     class_7.3-19       tools_4.1.1        bit_4.0.4         
 [9] memoise_2.0.0      RSQLite_2.2.8      lifecycle_1.0.1    tibble_3.1.5      
[13] lattice_0.20-45    pkgconfig_2.0.3    rlang_0.4.12       Matrix_1.3-4      
[17] igraph_1.2.7       DBI_1.1.1          rgdal_1.5-27       fastmap_1.1.0     
[21] e1071_1.7-9        terra_1.4-11       dplyr_1.0.7        hms_1.1.1         
[25] rgeos_0.5-8        generics_0.1.1     vctrs_0.3.8        bit64_4.0.5       
[29] classInt_0.4-3     grid_4.1.1         tidyselect_1.1.1   glue_1.4.2        
[33] R6_2.5.1           fansi_0.5.0        foreign_0.8-81     blob_1.2.2        
[37] purrr_0.3.4        magrittr_2.0.1     maptools_1.1-2     MASS_7.3-54       
[41] codetools_0.2-18   ellipsis_0.3.2     units_0.7-2        assertthat_0.2.1  
[45] utf8_1.2.2         KernSmooth_2.23-20 proxy_0.4-26       cachem_1.0.6      
[49] crayon_1.4.1       SSN_1.1.15