Skip to content

Instantly share code, notes, and snippets.

@h-a-graham

h-a-graham/weird_cutline_issue.md

Last active July 24, 2022 11:49
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save h-a-graham/7b6de5049dadc8b204e07fcc0b3f9300 to your computer and use it in GitHub Desktop.
Save h-a-graham/7b6de5049dadc8b204e07fcc0b3f9300 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
weird_cutline_issue.md 
library(vapour)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
#> WARNING: different compile-time and runtime versions for GEOS found:
#> Linked against: 3.10.2-CAPI-1.16.0 compiled against: 3.10.1-CAPI-1.16.0
#> It is probably a good idea to reinstall sf, and maybe rgeos and rgdal too
library(rstac)

mpc_dem <- function(aoi, src = c("cop-dem-glo-30", "alos-dem")) {
  stac_region <- sf::st_transform(aoi, 4326)  %>% 
    sf::st_bbox(aoi)
  
  it_obj <- rstac::stac("https://planetarycomputer.microsoft.com/api/stac/v1")|>
    rstac::stac_search(collections = src[1],
                       bbox = stac_region) |>
    rstac::get_request()
  
  src_list <- rstac::assets_list(it_obj)$path
  
  src_list[grep(".tif$", src_list)] 
}

f <- system.file("gpkg", "nc.gpkg", package = "sf")

src1 <- "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif"
src2 <- mpc_dem(read_sf(f)) 

info <- vapour::vapour_layer_info(f)

vapour::vapour_warp_raster(src1,
                           extent = info$extent,
                           dimension = c(10,10),
                           projection = info$projection$Wkt,
                           options=c("-cutline", f))
#> $Band1
#>   [1]        NaN        NaN        NaN 551.110229 342.492279 198.415161
#>   [7] 113.090706  69.613403  17.272562        NaN        NaN        NaN
#>  [13]        NaN 551.110229 342.492279 198.415161 113.090706  69.613403
#>  [19]  17.272562        NaN        NaN        NaN 821.525940 309.867645
#>  [25] 234.361374 175.136200  85.149147  43.532635  11.367164        NaN
#>  [31]        NaN        NaN 821.525940 309.867645 234.361374 175.136200
#>  [37]  85.149147  43.532635  11.367164        NaN 871.028870 953.713196
#>  [43] 401.690491 226.384598 194.823242 125.115105  51.705830  33.722717
#>  [49]   7.680205        NaN 871.028870 953.713196 401.690491 226.384598
#>  [55] 194.823242 125.115105  51.705830  33.722717   7.680205        NaN
#>  [61]        NaN        NaN        NaN        NaN        NaN  64.208275
#>  [67]  28.246149  18.623472  -4.679377        NaN        NaN        NaN
#>  [73]        NaN        NaN        NaN  64.208275  28.246149  18.623472
#>  [79]  -4.679377        NaN        NaN        NaN        NaN        NaN
#>  [85]        NaN        NaN  12.806995  -0.303200        NaN        NaN
#>  [91]        NaN        NaN        NaN        NaN        NaN        NaN
#>  [97]  12.806995  -0.303200        NaN        NaN

vapour::vapour_warp_raster(src2,
                           extent = info$extent,
                           dimension = c(10,10),
                           projection = info$projection$Wkt)
#> $Band1
#>   [1]  349.076111  411.445740 1123.195312 1039.526123  411.474243  180.004562
#>   [7]  149.584793   59.115612   16.175039    0.000000  322.887695  403.843292
#>  [13] 1117.813965  415.917084  274.266815  204.170700  119.703903   36.218086
#>  [19]   17.955286    0.000000  300.717285  367.891907 1012.624268  331.945343
#>  [25]  225.914169  211.283691   71.076454   34.534576    4.528758    0.000000
#>  [31]  540.184937  794.246643  467.625336  264.482544  215.672729  173.654694
#>  [37]   93.477135   40.137619   16.745296    1.578903  605.123352 1821.468506
#>  [43]  324.215210  249.290131  193.242279  101.210434   74.164200   34.413101
#>  [49]    8.329866    0.000000  931.524414  962.362305  260.718384  244.881897
#>  [55]  128.660110  109.053070   59.547546   34.601929    6.986561    0.000000
#>  [61]  656.597229  293.456512  220.223389  199.511612  162.050797   78.050270
#>  [67]   33.151157   23.180037   14.027452    0.000000  470.745453  246.343979
#>  [73]  220.954010  180.365143   99.621552   45.161091   24.663548    9.521521
#>  [79]    0.000000    0.000000  331.248199  233.221725  170.200378   96.720970
#>  [85]   99.985168   26.007788   20.302559    0.000000    0.000000    0.000000
#>  [91]  299.770966  165.619705  147.815781  124.719589   53.561546   18.991646
#>  [97]   18.455637    0.000000    0.000000    0.000000

vapour::vapour_warp_raster(src2,
                           extent = info$extent,
                           dimension = c(10,10),
                           projection = info$projection$Wkt,
                           options=c("-cutline", f))
#> Warning in warp_in_memory_gdal_cpp(x, source_WKT = source_projection, target_WKT
#> = projection, : GDAL Error 1: Cutline transformation failed
#> Error in warp_in_memory_gdal_cpp(x, source_WKT = source_projection, target_WKT = projection, : something went wrong!

Created on 2022-07-07 by the reprex package (v2.0.1)

@mdsumner
Copy link

mdsumner commented Jul 8, 2022

this works for me, bit slow because of the poor overview-scheme of the COP30 tiles, takes a few minutes to churn over every 27 tiles for a few pixels from the lowest overview - each tile is 1-degree so about 27 of them

prod(ceiling(diff(info$extent)[c(1, 3)]))
[1] 27
src2 <- "/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt"
vapour::vapour_warp_raster(src2,
                           extent = info$extent,
                           dimension = c(10,10),
                           projection = info$projection$Wkt,
                           options=c("-cutline", f))

$Band1
  [1]    0.000000    0.000000    0.000000 1037.989258  394.586060  178.896210  144.112396   55.927807   16.496422
 [10]    0.000000    0.000000    0.000000    0.000000  417.256409  269.472168  207.188553  115.525993   33.024612
 [19]   17.480848    0.000000    0.000000    0.000000 1012.866272  335.052490  207.444031  210.338486   60.022800
 [28]   31.684006   13.381020    0.000000    0.000000  751.996277  444.907837  272.380890  213.136490  182.453842
 [37]   92.997635   46.274082   18.117386    0.809912  623.721802 1785.558228  334.098480  242.366684  201.655991
 [46]   85.058609   71.958778   33.199013    8.000000    0.000000  890.873169  969.381226    0.000000    0.000000
 [55]  128.593048  108.265808   56.804798   38.164001   12.766655    0.000000    0.000000    0.000000    0.000000
 [64]    0.000000  168.818451   69.422043   32.613571   23.518702   10.010341    0.000000    0.000000    0.000000
 [73]    0.000000    0.000000    0.000000    0.000000   23.763950    8.916115    0.000000    0.000000    0.000000
 [82]    0.000000    0.000000    0.000000    0.000000    0.000000   16.891457    0.000000    0.000000    0.000000
 [91]    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000   20.906197    0.000000    0.000000
[100]    0.000000

@h-a-graham
Copy link
Author

I'd missed this, thankyou!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment