paul-carteron / happign

Happign allows you to use the APIs provided by the IGN (France) to download their public data.
https://paul-carteron.github.io/happign/
GNU General Public License v3.0
29 stars 2 forks source link

Add wms scan support #30

Open paul-carteron opened 1 day ago

paul-carteron commented 1 day ago
rm(list=ls())

library(sf)
#> Warning: package 'sf' was built under R version 4.3.3
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(terra)
#> Warning: package 'terra' was built under R version 4.3.3
#> terra 1.7.78

penmarch <- read_sf(system.file("extdata/penmarch.shp", package = "happign"))

get_scan <- function(x,
                     type = "scan25",
                     res = 10,
                     crs = 2154,
                     filename = tempfile(fileext = ".tif")){

   layer <- switch(type,
                   scan25 = "SCAN25TOUR_PYR-JPEG_WLD_WM",
                   scan100 = "SCAN100_PYR-JPEG_WLD_WM",
                   oaci = "SCANOACI_PYR-JPEG_WLD_WM",
                   carte_ign = "GEOGRAPHICALGRIDSYSTEMS.MAPS",
                   stop("layer should be one of 'scan25', 'scan100', 'oaci' or 'carte_ign'.",
                        call. = FALSE))

   wms_url <- sprintf("WMS:https://data.geopf.fr/private/wms-r?VERSION=1.3.0&LAYERS=%s&CRS=EPSG:4326&FORMAT=image/geotiff&apikey=ign_scan_ws",
                      layer)

   x <- x |>
      st_union() |>
      st_transform(crs)

   bb <- st_bbox(x)

   warp_options <- c(
      "-of", "GTIFF",
      "-te", bb$xmin, bb$ymin, bb$xmax, bb$ymax,
      "-te_srs", st_crs(crs)$srid,
      "-t_srs", st_crs(crs)$srid,
      "-tr", res, res,
      "-r", "bilinear",
      "-overwrite"
   )

   rast <- gdal_utils(
      util = "warp",
      source = wms_url,
      destination = filename,
      quiet = FALSE,
      options = c(warp_options)
   ) |> suppressWarnings()

   scan <- rast(filename)
   RGB(scan) <- c(1, 2, 3)
   names(scan) <- c("red", "green", "blue")

   return(scan)
}

scan25 <- get_scan(penmarch, type = "scan25", res = 1)
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
plotRGB(scan25)

Created on 2024-11-12 with reprex v2.1.1