Closed paul-carteron closed 12 months ago
gdal_translate with -tr and -projwin argument work but I need to compute width and height for each tile as before. Also I cannot reproject the raster.
rm(list=ls())
library(sf)
library(tmap)
library(terra)
# x <- mapedit::drawFeatures()
x <- read_sf(system.file("extdata/penmarch.shp", package = "happign"))
apikey = "ortho"
layer = "ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO"
res = 5
filename = tempfile(fileext = ".tif")
crs = 2154
overwrite = FALSE
version = "1.3.0"
styles = ""
interactive = FALSE
x <- st_make_valid(x) |>
st_transform(st_crs(crs))
bbox <- st_bbox(x)
url <- paste0(
"https://wxs.ign.fr/",
apikey,
"/geoportail/r/wms?",
"VERSION=1.3.0",
"&REQUEST=GetMap",
"&LAYERS=", layer,
"&CRS=EPSG:2154",
"&FORMAT=image/geotiff",
"&width=2048",
"&height=2048",
"&styles=",
"&bbox=", paste(bbox["xmin"], bbox["ymin"], bbox["xmax"], bbox["ymax"], sep = ",")
)
gdal_utils("translate",
source = url,
destination = "temp.tif",
quiet = FALSE,
options = c("-tr", 2, 2,
"-projwin", bbox$xmin, bbox$ymax, bbox$xmax, bbox$ymin))
rast <- rast("temp.tif")
tm_shape(rast)+
tm_rgb()
gdal warp can't be used because RGB is returned and also it's much longer
happign 0.1.9 work, need to go back to old method