resample misses some pixels

[carlos-arnillas]

Hello
I am using terra to reproject and combine several small rasters (+200). The code I created using terra worked mostly fine, but then I noticed gaps when merging some tiles (the versions of the installed packages that I'm using are listed below).
The code below downloads a tile from the Canadian DEM and reprojects and resample it (for simplicity, the resampling just estimates the total number of pixel). The original raster is in ~80x80m and the resampling transforms it into 1x1km in a UTM coordinates. After resampling, the code creates a figure showing the mismatch (pasted below). In the figure, the dashed red lines are the left and top boundary of the coarse raster. Notice how the Northern section of the 80x80m raster is not transformed, staying as NA values (hence transparent) in the coarse resolution raster.
(a similar problem occurs in other raster files, but the missing section could be the bottom row, instead of the top one, so this is just an example).
Thank you for your help.

Carlos Alberto

library(terra)
library(sf)

# Loading the data
download.file(paste0("https://ftp.maps.canada.ca/pub/nrcan_rncan/archive/elevation/geobase_cded_dnec/250k_dem/057/057c.zip"), 
              paste0("x057c.zip"))
unzip("x057c.zip")
crs_planar <- "EPSG:32615"  # UTM 15N

# Reading the info
r1 <- rast("057c_0100_deme.dem")
# reprojecting to UTM
r2 <- project(r1, crs_planar)
# building the bbox
bb <- st_bbox(r2)
bb[c("xmin", "ymin")] <- floor(bb[c("xmin", "ymin")] / 1000)*1000 - 500
bb[c("xmax", "ymax")] <- ceiling(bb[c("xmax", "ymax")] / 1000)*1000 + 500

# creating the reference raster
gr <- rast(resolution = 1000, 
           xmin=bb["xmin"], xmax=bb["xmax"], ymin=bb["ymin"], ymax=bb["ymax"])
# and resampling
rn <- resample(!is.na(r2), gr, "sum")

# final plots
plot(r2, xlim=c(455000, 470000), ylim=c(7760000,7770000))
abline(a=bb["ymax"], b=0, col="red", lty="dashed")
abline(v=bb["xmin"], col="red", lty="dashed")
plot(rn, xlim=c(455000, 470000), ylim=c(7760000,7770000), add=T, legend=F)
# notice the mismatch in the areas at the top: at least one more 1x1k pixel should have been estimated by resample
# so that in the final plot the fine resolution raster should not be visible at all.

This is the current installation (they are the latest available when I try installing them in my Ubuntu 20.04 computer):

> terra::gdal(lib="all")
    gdal     proj     geos 
 "3.4.3"  "8.2.0" "3.10.2" 
> packageVersion("terra")
[1] ‘1.7.83’

[rhijmans]

This works as expected for me with the CRAN version of "terra" on windows.

library(terra)
#terra 1.7.83

if (!file.exists("057c_0100_deme.dem")) {
	url <- "https://ftp.maps.canada.ca/pub/nrcan_rncan/archive/elevation/geobase_cded_dnec/250k_dem/057/057c.zip"
	download.file(url, basename(url), mode="wb") 
	unzip(basename(url))
}

r1 <- rast("057c_0100_deme.dem")
r2 <- project(r1, "EPSG:32615") # UTM 15N
e <- round(ext(r2), -3) + 1500
gr <- rast(e, res=1000, crs=crs(r2))
rn <- resample(!is.na(r2), gr, "sum")

plot(r2, xlim=c(455000, 470000), ylim=c(7760000,7770000))
abline(a=e$ymax, b=0, col="red", lty="dashed")
abline(v=e$xmin, col="red", lty="dashed")
plot(rn, xlim=c(455000, 470000), ylim=c(7760000,7770000), add=T, legend=F)

gdal(lib="all")
#    gdal     proj     geos 
# "3.9.3"  "9.4.1" "3.12.2" 

[rhijmans]

It also works as expected with the development version on Ubuntu 24.04.1 LTS with

 gdal(lib="all")
 #   gdal     proj     geos
 #"3.9.3"  "9.4.1" "3.12.2"

So perhaps it is a matter of the GDAL version.

[carlos-arnillas]

thank you I will try updating all packages.

…

On Mon, Nov 25, 2024 at 9:40 PM Robert Hijmans ***@***.***> wrote: It also works as expected with the development version on Ubuntu 24.04.1 LTS with gdal(lib="all") # gdal proj geos #"3.9.3" "9.4.1" "3.12.2" So perhaps it is a matter of the GDAL version. — Reply to this email directly, view it on GitHub <#1651 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGGQELPAQNMRN7OGVQR7TXL2CPNQNAVCNFSM6AAAAABSOXVMCCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIOJZGUZDQMZSG4> . You are receiving this because you authored the thread.Message ID: ***@***.***>