Crop with mask=TRUE reduces all values to 1, and saving to memory results in : Error: [crop] file exists. You can use 'overwrite=TRUE' to overwrite it

[Rapsodia86]

I do have a big raster (72439 x 88523 px). When trying to crop with mask=TRUE, all values are turned to 1. It is like, instead of the cropped raster, I was getting the mask (rasterized polygon). If I split the process into two steps (first crop, then mask), everything is ok. If I save it to a file, instead of to memory, the problem persists. With a smaller shapefile everything is ok. However, when wanted to rerun the crop() with mask=TRUE, I got an error: [crop] file exists. You can use 'overwrite=TRUE' to overwrite it

All steps below:

> library(terra)
terra 1.7.74
Warning message:
package ‘terra’ was built under R version 4.2.3 
> eucrop18 <- rast("MyDir/EUCROPMAP_2018_crop_3844_verR.tif")*1
> eucrop18                                
class       : SpatRaster 
dimensions  : 72439, 88523, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 57543.57, 942773.6, 140629.4, 865019.4  (xmin, xmax, ymin, ymax)
coord. ref. : Pulkovo 1942(58) / Stereo70 (EPSG:3844) 
source      : spat_342c447c6553_13356.tif 
varname     : EUCROPMAP_2018_crop_3844_verR 
name        : EUCROPMAP_2018_crop_3844 
min value   :                      100 
max value   :                      600 
> plot(eucrop18)

> shp_big <- vect("MyDir/Unitate_administrativa_judet_3844_all.shp")
> shp_big
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 42, 17  (geometries, attributes)
 extent      : 134104, 870647.1, 235538.6, 753220.1  (xmin, xmax, ymin, ymax)
 source      : Unitate_administrativa_judet_3844_all.shp
 coord. ref. : Pulkovo 1942(58) / Stereo70 (EPSG:3844) 
 names       : OBJECTID    Id IFCID localId   namespace      versionId
 type        :    <int> <int> <int>   <chr>       <chr>          <chr>
 values      :        1  3186  3186    1.29 RO.Ancpi.AU 20230116112626
                      2  3187  3187    1.38 RO.Ancpi.AU 20230116112626
                      3  3188  3188    1.47 RO.Ancpi.AU 20230116112627
 country natLevel natLevName natCode (and 7 more)
   <chr>    <chr>      <chr>   <chr>             
      RO 2ndOrder      Judet      29             
      RO 2ndOrder      Judet      38             
      RO 2ndOrder      Judet      47   
> plot(shp_big,add=TRUE)

crop_big <- crop(eucrop18,shp_big,mask=TRUE,touches=FALSE) 
plot(crop_big)

Two-step approach:

> crop_big <- crop(eucrop18,shp_big,touches=FALSE) 
> plot(crop_big)

crop_big_mask <- mask(crop_big,shp_big,touches=FALSE) 
plot(crop_big_mask)

Small shapefile scenario:

shp_small <- vect("MyDir/Unitate_administrativa_judet_3844_sel5_WEST.shp")
> shp_small
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 5, 17  (geometries, attributes)
 extent      : 134104, 412619.9, 416432.1, 737281.1  (xmin, xmax, ymin, ymax)
 source      : Unitate_administrativa_judet_3844_sel5_WEST.shp
 coord. ref. : Pulkovo 1942(58) / Stereo70 (EPSG:3844) 
 names       : OBJECTID    Id IFCID localId   namespace      versionId
 type        :    <int> <int> <int>   <chr>       <chr>          <chr>
 values      :        1  3186  3186    1.29 RO.Ancpi.AU 20230116112626
                      4  3189  3189    1.56 RO.Ancpi.AU 20230116112627
                     18  3203  3203   1.305 RO.Ancpi.AU 20230116112628
 country natLevel natLevName natCode (and 7 more)
   <chr>    <chr>      <chr>   <chr>             
      RO 2ndOrder      Judet      29             
      RO 2ndOrder      Judet      56             
      RO 2ndOrder      Judet     305             
> plot(eucrop18)
> plot(shp_small,add=TRUE)

> crop_small <- crop(eucrop18,shp_small,mask=TRUE,touches=FALSE) 
> plot(crop_small)

Finally, when wanted to repeat the process with mask=TRUE, I got the error.

> crop_big2 <- crop(eucrop18,shp_big,mask=TRUE,touches=FALSE) 
Error: [crop] file exists. You can use 'overwrite=TRUE' to overwrite it

Session info:

sessionInfo() R version 4.2.2 (2022-10-31 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.utf8 [2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

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

other attached packages: [1] terra_1.7-74

loaded via a namespace (and not attached): [1] Rcpp_1.0.12 codetools_0.2-18 fansi_1.0.6 utf8_1.2.4
[5] dplyr_1.1.4 parallelly_1.37.1 plyr_1.8.9 R6_2.5.1
[9] lifecycle_1.0.4 magrittr_2.0.3 pillar_1.9.0 stringi_1.8.3
[13] rlang_1.1.3 cli_3.6.2 reshape2_1.4.4 rstudioapi_0.15.0 [17] vctrs_0.6.5 generics_0.1.3 tools_4.2.2 stringr_1.5.1
[21] glue_1.7.0 parallel_4.2.2 compiler_4.2.2 pkgconfig_2.0.3
[25] tidyselect_1.2.1 tibble_3.2.1

[kuntalnr]

Hello. I am also facing the same issue. When I crop and mask using terra::crop(spat_raster, shp, snap = "out", mask = T), I have to add overwrite = T and the output files are all converted to 1. The terra version I am using is 1.7.71

[Rapsodia86]

That raster with value of 1 looks like it is the rasterized polygon that is used as a mask.

[kuntalnr]

Hello. Just checking if this issue got resolved?

[kadyb]

As a workaround, you can try using gdalwarp with, for example, sf::gdal_utils():

library("sf")
library("terra")

tmp = rast(ncols = 100, nrows = 100, xmin = 5.74414, xmax = 6.528252,
           ymin = 49.44781, ymax = 50.18162, crs = "EPSG:4326")
init(tmp, fun = 1, filename = "test.tif")
f = system.file("ex/lux.shp", package = "terra")

sf::gdal_utils(util = "warp", source = "test.tif", destination = "output.tif",
               options = c("-cutline", f, "-crop_to_cutline"))
rast("output.tif")

[Rapsodia86]

Rather just do this in two steps (two functions): first crop, then mask the results. At least you will get the same result if the function was working correctly. Using another package may result in a diifferent clipping/masking approach (e.g. you may end up with different pixels along the border) if there is no control about this.

[ailich]

Also running into this issue after updating terra

[Rapsodia86]

Hi @rhijmans ! Could you take a look to check what is going on with crop and mask functions when you have a moment? Appreciate very much! Monika

[rhijmans]

Thank you for your patience. I believe this now also works for SpatRasters that are written to file.

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
v <- vect(system.file("ex/lux.shp", package="terra"))
m <- crop(r, v[9:12,], mask=TRUE, todisk=TRUE)
plot(m)

[Rapsodia86]

Thank you very much!

[Rimagination]

I have recently encountered this problem.

[Rapsodia86]

I have recently encountered this problem.

Are you using the most recent version of the package?