search

opengeos / whiteboxR

WhiteboxTools R Frontend
https://whiteboxR.gishub.org
Other
172 stars 30 forks source link

wbt_multiscale_topographic_position_image() throws error despite input rasters having same dimensions and extent #118

Closed evhersh closed 11 months ago

evhersh commented 11 months ago

Thanks for the great package! I ran into an issue while trying to make a multiscale topographic position image. I first made 3 DEVmax rasters using wbt_max_elevation_deviation() with min/max scales of 5-20, 20-200, and 200-2000 all based on the same input DEM. After trying to run the MTP function, I receive this error: 

> wbt_multiscale_topographic_position_image(local = here("AP_DEM", "devmax5-20.tif"),
+                                           meso = here("AP_DEM", "devmax20-200.tif"),
+                                           broad = here("AP_DEM", "devmax200-2000.tif"),
+                                           output = here("AP_DEM", "MTP_5-2000.tif"))

Error running WhiteboxTools (MultiscaleTopographicPositionImage)
  whitebox.exe_path: "C:\Users\evanh\AppData\Roaming/R/data/R/whitebox/WBT/whitebox_tools.exe"; File exists? TRUE
  Arguments: --run=MultiscaleTopographicPositionImage  --local="C:/GitHub/Galore_PEM/AP_DEM/devmax5-20.tif" --meso="C:/GitHub/Galore_PEM/AP_DEM/devmax20-200.tif" --broad="C:/GitHub/Galore_PEM/AP_DEM/devmax200-2000.tif" --output="C:/GitHub/Galore_PEM/AP_DEM/MTP_5-2000.tif" --lightness=1.2 -v

System command had status 101
*************************************************
* Welcome to MultiscaleTopographicPositionImage *
* Powered by WhiteboxTools                      *
* www.whiteboxgeo.com                           *
*************************************************
Reading broad-scale DEV data...
Reading meso-scale DEV data...
Reading local-scale DEV data...
thread 'main' panicked at 'The input files must have the same number of rows and columns and spatial extent.', whitebox-tools-app\src\main.rs:72:21
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
multiscale_topographic_position_image - Elapsed Time: NA [did not run]

I checked the raster dimensions and extents using the terra package, and they all appear to have identical values:

> rast(here("AP_DEM", "devmax5-20.tif") )
class       : SpatRaster 
dimensions  : 13396, 15663, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 303100, 459730, 6248960, 6382920  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : devmax5-20.tif 
name        : devmax5-20 

> rast(here("AP_DEM", "devmax20-200.tif"))
class       : SpatRaster 
dimensions  : 13396, 15663, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 303100, 459730, 6248960, 6382920  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : devmax20-200.tif 
name        : devmax20-200 

> rast(here("AP_DEM", "devmax200-2000.tif"))
class       : SpatRaster 
dimensions  : 13396, 15663, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 303100, 459730, 6248960, 6382920  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : devmax200-2000.tif 
name        : devmax200-2000

Perhaps there is something I'm missing? I am using package version 2.3.3 in RStudio on a machine running Windows 11. All of the other functions have been working great up until now.

Thanks for your help!

Evan

evhersh commented 11 months ago

Hello again,

I just wanted to follow up on this with some more information. I've tested this on my Macbook using the same exact code and input files and it is working as expected, so this must be related to the Windows version. I tested it with the example DEM that came with the package, and confirmed that the issue still exists on my Windows machine (so it is not related to peculiarities with my input DEM). 

library(terra)
library(here)
library(whitebox)

wbt_init()

dem <- file.path(system.file("extdata", package = "whitebox"), "DEM.tif")
plot(rast(dem))

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 5,
                            max_scale = 20,
                            out_mag = here("test_devmax5-20.tif"),
                            out_scale= here("test_devscale5-20.tif"))

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 20,
                            max_scale = 200,
                            out_mag = here("test_devmax20-200.tif"),
                            out_scale= here("test_devscale20-200.tif"))

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 200,
                            max_scale = 2000,
                            out_mag = here("test_devmax200-2000.tif"),
                            out_scale= here("test_devscale200-2000.tif"))

wbt_multiscale_topographic_position_image(local = here("test_devmax5-20.tif"),
                                          meso = here("test_devmax20-200.tif"),
                                          broad = here("test_devmax200-2000.tif"),
                                          output = here("test_MTP_5-2000.tif"))

Error running WhiteboxTools (MultiscaleTopographicPositionImage)
  whitebox.exe_path: "C:\Users\evanh\AppData\Roaming/R/data/R/whitebox/WBT/whitebox_tools.exe"; File exists? TRUE
  Arguments: --run=MultiscaleTopographicPositionImage  --local="C:/GitHub/Lawyers_PEM/test_devmax5-20.tif" --meso="C:/GitHub/Lawyers_PEM/test_devmax20-200.tif" --broad="C:/GitHub/Lawyers_PEM/test_devmax200-2000.tif" --output="C:/GitHub/Lawyers_PEM/test_MTP_5-2000.tif" --lightness=1.2 -v

System command had status 101
*************************************************
* Welcome to MultiscaleTopographicPositionImage *
* Powered by WhiteboxTools                      *
* www.whiteboxgeo.com                           *
*************************************************
Reading broad-scale DEV data...
Reading meso-scale DEV data...
Reading local-scale DEV data...
thread 'main' panicked at 'The input files must have the same number of rows and columns and spatial extent.', whitebox-tools-app\src\main.rs:72:21
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
multiscale_topographic_position_image - Elapsed Time: NA [did not run]
brownag commented 11 months ago

Hi @evhersh

Thanks for reporting and providing the reproducible example. I can confirm the 'The input files must have the same number of rows and columns and spatial extent.' error occurs on my Linux machine.

If I supply the optional hillshade argument, like below, everything "works", though there are warnings issued due to the small dimensions in the demo raster v.s. your selected min/max scale parameters:

library(terra)
#> terra 1.7.49
library(here)
library(whitebox)

wbt_init(verbose = TRUE)

dem <- file.path(system.file("extdata", package = "whitebox"), "DEM.tif")

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 5,
                            max_scale = 20,
                            out_mag = here("test_devmax5-20.tif"),
                            out_scale= here("test_devscale5-20.tif"))
#> max_elevation_deviation - Elapsed Time (excluding I/O): 0.19s

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 20,
                            max_scale = 200,
                            out_mag = here("test_devmax20-200.tif"),
                            out_scale= here("test_devscale20-200.tif"))
#> Warning: The number of steps resulted in filter sizes that 
#> max_elevation_deviation - Elapsed Time (excluding I/O): 0.26s

wbt_max_elevation_deviation(dem=dem,
                            min_scale = 200,
                            max_scale = 2000,
                            out_mag = here("test_devmax200-2000.tif"),
                            out_scale= here("test_devscale200-2000.tif"))
#> Warning: The number of steps resulted in filter sizes that 
#> max_elevation_deviation - Elapsed Time (excluding I/O): 0.1s

terra::describe(here("test_devmax5-20.tif"))
#>  [1] "Driver: GTiff/GeoTIFF"                                                                                                                                                                                                                                                          
#>  [2] "Files: /tmp/RtmpbO0xoH/reprex-71d576087f3f-cummy-krill/test_devmax5-20.tif"                                                                                                                                                                                                     
#>  [3] "Size is 237, 188"                                                                                                                                                                                                                                                               
#>  [4] "Coordinate System is:"                                                                                                                                                                                                                                                          
#>  [5] "PROJCRS[\"NAD83 / UTM zone 18N\","                                                                                                                                                                                                                                              
#>  [6] "    BASEGEOGCRS[\"NAD83\","                                                                                                                                                                                                                                                     
#>  [7] "        DATUM[\"North American Datum 1983\","                                                                                                                                                                                                                                   
#>  [8] "            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,"                                                                                                                                                                                                                      
#>  [9] "                LENGTHUNIT[\"metre\",1]]],"                                                                                                                                                                                                                                     
#> [10] "        PRIMEM[\"Greenwich\",0,"                                                                                                                                                                                                                                                
#> [11] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                                                                                                                                                                                                                         
#> [12] "        ID[\"EPSG\",4269]],"                                                                                                                                                                                                                                                    
#> [13] "    CONVERSION[\"UTM zone 18N\","                                                                                                                                                                                                                                               
#> [14] "        METHOD[\"Transverse Mercator\","                                                                                                                                                                                                                                        
#> [15] "            ID[\"EPSG\",9807]],"                                                                                                                                                                                                                                                
#> [16] "        PARAMETER[\"Latitude of natural origin\",0,"                                                                                                                                                                                                                            
#> [17] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                                                                                                                                                                                                                          
#> [18] "            ID[\"EPSG\",8801]],"                                                                                                                                                                                                                                                
#> [19] "        PARAMETER[\"Longitude of natural origin\",-75,"                                                                                                                                                                                                                         
#> [20] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                                                                                                                                                                                                                          
#> [21] "            ID[\"EPSG\",8802]],"                                                                                                                                                                                                                                                
#> [22] "        PARAMETER[\"Scale factor at natural origin\",0.9996,"                                                                                                                                                                                                                   
#> [23] "            SCALEUNIT[\"unity\",1],"                                                                                                                                                                                                                                            
#> [24] "            ID[\"EPSG\",8805]],"                                                                                                                                                                                                                                                
#> [25] "        PARAMETER[\"False easting\",500000,"                                                                                                                                                                                                                                    
#> [26] "            LENGTHUNIT[\"metre\",1],"                                                                                                                                                                                                                                           
#> [27] "            ID[\"EPSG\",8806]],"                                                                                                                                                                                                                                                
#> [28] "        PARAMETER[\"False northing\",0,"                                                                                                                                                                                                                                        
#> [29] "            LENGTHUNIT[\"metre\",1],"                                                                                                                                                                                                                                           
#> [30] "            ID[\"EPSG\",8807]]],"                                                                                                                                                                                                                                               
#> [31] "    CS[Cartesian,2],"                                                                                                                                                                                                                                                           
#> [32] "        AXIS[\"(E)\",east,"                                                                                                                                                                                                                                                     
#> [33] "            ORDER[1],"                                                                                                                                                                                                                                                          
#> [34] "            LENGTHUNIT[\"metre\",1]],"                                                                                                                                                                                                                                          
#> [35] "        AXIS[\"(N)\",north,"                                                                                                                                                                                                                                                    
#> [36] "            ORDER[2],"                                                                                                                                                                                                                                                          
#> [37] "            LENGTHUNIT[\"metre\",1]],"                                                                                                                                                                                                                                          
#> [38] "    USAGE["                                                                                                                                                                                                                                                                     
#> [39] "        SCOPE[\"Engineering survey, topographic mapping.\"],"                                                                                                                                                                                                                   
#> [40] "        AREA[\"North America - between 78°W and 72°W - onshore and offshore. Canada - Nunavut; Ontario; Quebec. United States (USA) - Connecticut; Delaware; Maryland; Massachusetts; New Hampshire; New Jersey; New York; North Carolina; Pennsylvania; Virginia; Vermont.\"],"
#> [41] "        BBOX[28.28,-78,84,-72]],"                                                                                                                                                                                                                                               
#> [42] "    ID[\"EPSG\",26918]]"                                                                                                                                                                                                                                                        
#> [43] "Data axis to CRS axis mapping: 1,2"                                                                                                                                                                                                                                             
#> [44] "Origin = (664692.000000000000000,4895824.000000000000000)"                                                                                                                                                                                                                      
#> [45] "Pixel Size = (90.000000000000000,-90.000000000000000)"                                                                                                                                                                                                                          
#> [46] "Metadata:"                                                                                                                                                                                                                                                                      
#> [47] "  AREA_OR_POINT=Area"                                                                                                                                                                                                                                                           
#> [48] "  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)"                                                                                                                                                                                                                                       
#> [49] "  TIFFTAG_SOFTWARE=WhiteboxTools"                                                                                                                                                                                                                                               
#> [50] "  TIFFTAG_XRESOLUTION=72"                                                                                                                                                                                                                                                       
#> [51] "  TIFFTAG_YRESOLUTION=72"                                                                                                                                                                                                                                                       
#> [52] "Image Structure Metadata:"                                                                                                                                                                                                                                                      
#> [53] "  INTERLEAVE=BAND"                                                                                                                                                                                                                                                              
#> [54] "Corner Coordinates:"                                                                                                                                                                                                                                                            
#> [55] "Upper Left  (  664692.000, 4895824.000) ( 72d56'20.40\"W, 44d11'49.36\"N)"                                                                                                                                                                                                      
#> [56] "Lower Left  (  664692.000, 4878904.000) ( 72d56'39.44\"W, 44d 2'41.31\"N)"                                                                                                                                                                                                      
#> [57] "Upper Right (  686022.000, 4895824.000) ( 72d40'20.20\"W, 44d11'30.91\"N)"                                                                                                                                                                                                      
#> [58] "Lower Right (  686022.000, 4878904.000) ( 72d40'41.69\"W, 44d 2'22.95\"N)"                                                                                                                                                                                                      
#> [59] "Center      (  675357.000, 4887364.000) ( 72d48'30.42\"W, 44d 7' 6.41\"N)"                                                                                                                                                                                                      
#> [60] "Band 1 Block=237x1 Type=Float32, ColorInterp=Gray"                                                                                                                                                                                                                              
#> [61] "  NoData Value=-32768"
# terra::describe(here("test_devmax20-200.tif"))
# terra::describe(here("test_devmax200-2000.tif"))

wbt_multidirectional_hillshade(dem, here("hillshade.tif"))
#> multidirectional_hillshade - Elapsed Time (excluding I/O): 0.5s

wbt_multiscale_topographic_position_image(local = here("test_devmax5-20.tif"),
                                          meso = here("test_devmax20-200.tif"),
                                          broad = here("test_devmax200-2000.tif"), 
                                          hillshade = here("hillshade.tif"),
                                          output = here("test_MTP_5-2000.tif"))
#> multiscale_topographic_position_image - Elapsed Time (excluding I/O): 0.1s
plotRGB(rast(here("test_MTP_5-2000.tif")))

It appears that this is not an issue with the number of rows/columns/spatial extent in the input rasters, as you point out they are identical, but rather some logic in lieu of the optional hillshade input that does not quite work right. Since the error does not occur when hillshade is specified, all other inputs equal, I estimate the issue is occurring around here: https://github.com/jblindsay/whitebox-tools/blob/a0746ecb0340eb7be6401c555f1629ebfa5d2de6/whitebox-tools-app/src/tools/terrain_analysis/multiscale_topographic_position_image.rs#L349-L356

As far as I can tell this is not an issue with the R package frontend, but rather an internal problem with the Multiscale Topographic Position Image tool. We can't fix the tool from here, so please report the over on https://github.com/jblindsay/whitebox-tools/issues