Add gdal_footprint support to gdal_utils

[goergen95]

Hi,
as discussed in #2304 this PR adds support for the gdal_footprint utility requiring GDAL >= 3.8.0.
I am very new to C++, so please double-check if everything is ok. I used the script below to check the functionality:

# remotes::install_github("https://github.com/goergen95/sf", ref = "gdalfootprint")
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.8.1, PROJ 9.1.1; sf_use_s2() is TRUE
tif <- system.file("tif/geomatrix.tif", package = "sf")
t1 <- tempfile(fileext = ".gpkg")
t2 <- tempfile(fileext = ".json")
t3 <- tempfile(fileext = ".gpkg")
t4 <- tempfile(fileext = ".gpkg")
t5 <- tempfile(fileext = ".gpkg")
t6 <- tempfile(fileext = ".gpkg")

# GPKG output
gdal_utils("gdalfootprint", tif, t1)
(s1 <- st_read(t1))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf77437ef6.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#>                             geom
#> 1 MULTIPOLYGON (((1841002 114...
# GeoJSON output
gdal_utils("gdalfootprint", tif, t2, options = c("-of", "GeoJSON"))
(s2 <- st_read(t2))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf1a8c155b.json' using driver `GeoJSON'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#>                         geometry
#> 1 MULTIPOLYGON (((1841002 114...
# gdal_footprint options
gdal_utils("gdalfootprint", tif, t3, options = c("-t_srs", "EPSG:4326"))
(s3 <- st_read(t3))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf33b7643b.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -104.8474 ymin: 10.11931 xmax: -104.8463 ymax: 10.12043
#> Geodetic CRS:  WGS 84
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -104.8474 ymin: 10.11931 xmax: -104.8463 ymax: 10.12043
#> Geodetic CRS:  WGS 84
#>                             geom
#> 1 MULTIPOLYGON (((-104.8465 1...
# opening options
gdal_utils("gdalfootprint", tif, t4, options = c("-oo", "NUM_THREADS=1"))
(s4 <- st_read(t4))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf1086b216.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#>                             geom
#> 1 MULTIPOLYGON (((1841002 114...
# data set creation options
gdal_utils("gdalfootprint", tif, t5, options = c("-dsco", "CRS_WKT_EXTENSION=Yes"))
(s5 <- st_read(t5))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf1d7b91f7.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1840902 ymin: 1143873 xmax: 1841032 ymax: 1144003
#> Projected CRS: WGS 84 / UTM zone 11N
#>                             geom
#> 1 MULTIPOLYGON (((1841002 114...
# remote file with overview option
tif2 <- paste0("/vsicurl/https://deafrica-input-datasets.s3.af-south-1.amazonaws.com/", 
               "cci_landcover/1992/cci_landcover_1992_v2.0.7cds.tif")
gdal_utils("gdalfootprint", tif2, t6, options = c("-ovr", "4"))
(s6 <- st_read(t6))
#> Reading layer `footprint' from data source `/tmp/RtmpS1pM0Y/file27bf6fbdb9e1.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -26.36111 ymin: -47.96944 xmax: 64.5 ymax: 38.35
#> Geodetic CRS:  WGS 84
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -26.36111 ymin: -47.96944 xmax: 64.5 ymax: 38.35
#> Geodetic CRS:  WGS 84
#>                             geom
#> 1 MULTIPOLYGON (((-26.36111 3...

^{Created on 2024-01-01 with reprex v2.0.2}

[edzer]

This needs to be an #else, and the #endif should go after line 284.

[goergen95]

Thanks! I think I now use #if properly.

[edzer]

I'd suggest to call the util footprint, rather than gdalfootprint, in line with the others. Also above in the doc section.

[goergen95]

Done!

[edzer]

Great, thanks!

[edzer]

I added a few more commits to clean up.

[goergen95]

Cool, thanks for the help. I'll close the issue then.