the in-memory raster warper

mdsumner commented 4 years ago

[x] allow empty source_WKT
[ ] ensure empty and missing source proj is not a showstopper
[ ] allow projstring input
[ ] allow target proj to be empty, so it's just a RasterIO

## target coordinate system
  tas_wkt <- "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"unknown\",\n        METHOD[\"Lambert Azimuthal Equal Area\",\n            ID[\"EPSG\",9820]],\n        PARAMETER[\"Latitude of natural origin\",-42,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",147,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
prj <- "+proj=laea +lon_0=147 +lat_0=-42 +datum=WGS84"

## source coordinate system (GDAL often can't get one from a source, and needs a proxy cache for this to work by user-input)
ll_wkt <- "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

## arbitrary scaling factor
fac <- 1
dm <- as.integer(c(512, 298) * 3 * fac)  ## hardcoded for this example
gt <- c(-637239.4, 5030.0/fac, 0.0, 261208.7, 0.0, -7760.0/fac) ## geotransform version of extent+dim

gt_dim_2ex <- function(x, dim) {
  xx <- c(x[1], x[1] + dim[1] * x[2]) - x[2]/2
  yy <- c(x[4] + dim[2] * x[6], x[4]) - x[6]/2
  raster::extent(c(xx, yy))
}

ex <- gt_dim_2ex(gt, dm)  ## extent from geotransform
f <- system.file("extdata", "sst.tif", package = "vapour")
f <- raadtools::topofile("etopo2")  ## local copy of a longlat raster of the globe
## warp the source grid etopo2 to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(f, source_WKT = ll_wkt,target_geotransform = gt, target_dim = dm, target_WKT = tas_wkt, band = 1)
#> setting projection

## turn that into something we can visualize                             
v <- vals[[1]]
#v[v < 250] <- NA
library(raster)
#> Loading required package: sp
data("wrld_simpl", package = "maptools")
plot(setValues(raster(ex, nrows = dm[2], ncols = dm[1]), v))
plot(spTransform(wrld_simpl, prj), add = TRUE)

^{Created on 2020-05-28 by the reprex package (v0.3.0)}

mdsumner commented 4 years ago

consider that a warper interface could act as a RasterIO replace also, might be smarter than using the vapour_read_raster_cpp

mdsumner commented 4 years ago

more generic example (kinda), need that geotransform / extent-dim translator from affinity

## target coordinate system
prj <- "+proj=laea +lon_0=0 +lat_0=-90 +datum=WGS84"
wkt <- sf::st_crs(prj)$wkt
library(measoshapes)
p_mea <- subset(measo_regions05, !grepl("T$", name))

## source coordinate system (GDAL often can't get one from a source, and needs a proxy cache for this to work by user-input)
ll_wkt <- "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

dm <- c(1024, 1024)
target <- raster::raster(p_mea, nrows = dm[1], ncols = dm[2], crs = prj)
target[] <- 0
tf <- tempfile(fileext = ".tif")
raster::writeRaster(target, tf)
gt <- vapour::vapour_raster_info(tf)$geotransform

gt_dim_2ex <- function(x, dim) {
  xx <- c(x[1], x[1] + dim[1] * x[2]) - x[2]/2
  yy <- c(x[4] + dim[2] * x[6], x[4]) - x[6]/2
  raster::extent(c(xx, yy))
}

ex <- gt_dim_2ex(gt, dm)  ## extent from geotransform

f <- raadtools::topofile("gebco_19")  ## local copy of a longlat raster of the globe
## warp the source grid to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(f, source_WKT = ll_wkt,target_geotransform = gt, target_dim = dm, target_WKT = wkt, band = 1)

nn <- 24
image(setValues(target, vals[[1]]), 
      col = grey(seq(0, 1, length.out = nn -1)),
      asp = 1,
      breaks = quantile(vals[[1]], seq(0, 1, length.out = nn)))

mdsumner commented 3 years ago

Bit simpler now using functions in {affinity}

## target coordinate system
prj <- "+proj=laea +lon_0=0 +lat_0=-90 +datum=WGS84"
wkt <- sf::st_crs(prj)$wkt
library(measoshapes)
p_mea <- subset(measo_regions05, !grepl("T$", name))

## source coordinate system (GDAL often can't get one from a source, and needs a proxy cache for this to work by user-input)
ll_wkt <- "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"
dm <- c(1024, 1024)

target <- raster::raster(p_mea, nrows = dm[1], ncols = dm[2], crs = prj)
target[] <- 0
gt <- affinity:::raster_to_gt(target)
## extent from geotransform
ex <- affinity::gt_dim_to_extent(gt, dm)

f <- raadtools::topofile("gebco_19")  ## local copy of a longlat raster of the globe
## warp the source grid to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(f, source_WKT = ll_wkt,
                                            target_geotransform = gt, 
                                            target_dim = dm, 
                                            target_WKT = wkt, 
                                            band = 1)

nn <- 24
image(setValues(target, vals[[1]]), 
      col = grey(seq(0, 1, length.out = nn -1)),
      asp = 1,
      breaks = quantile(vals[[1]], seq(0, 1, length.out = nn)))

mdsumner commented 3 years ago

This example a little more relevant, note there's no need to set source_WKT:

## target coordinate system
prj <- "+proj=laea +lon_0=0 +lat_0=-90 +datum=WGS84"
wkt <- sf::st_crs(prj)$wkt

## source coordinate system (GDAL often can't get one from a source, and needs a proxy cache for this to work by user-input)
src_wkt <- ""
dm <- c(1024, 1024) 

target <- raster::raster(extent(-3e6, 3.2e6, -2.3e6, 2.3e6)/2, nrows = dm[1], ncols = dm[2], crs = prj)
target[] <- 0
gt <- affinity:::raster_to_gt(target)
## extent from geotransform
ex <- affinity::gt_dim_to_extent(gt, dm)

library(raadtools)
f <- raadfiles::rema_200m_files()$fullname[1]

## warp the source grid to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(f, source_WKT = src_wkt,
                                            target_geotransform = gt, 
                                            target_dim = dm, 
                                            target_WKT = wkt, 
                                            band = 1)

nn <- 24
set_na <- function(x) {
  x[x < -9000] <- NA
  x
}
vv <- set_na(vals[[1]])
image(setValues(target, vv), 
      col = grey(seq(0, 1, length.out = nn -1)),
      asp = 1,
      breaks = quantile(vv[vv > 0], seq(0, 1, length.out = nn), na.rm = TRUE))

mdsumner commented 3 years ago

a bit more interesting with online server config, but not working

https://gist.github.com/mdsumner/91b2dfe5c7acba9b3aa8fb30d01b8bad

mdsumner commented 3 years ago

WAAAAY better now

  prj <- "+proj=laea +lon_0=0 +lat_0=-90 +datum=WGS84"
wkt <- sf::st_crs(prj)$wkt
library(measoshapes)
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
p_mea <- subset(measo_regions05, !grepl("T$", name))

## source coordinate system (GDAL often can't get one from a source, and needs a proxy cache for this to work by user-input)
ll_wkt <- "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

dm <- c(1024, 1024)
target <- raster::raster(p_mea, nrows = dm[1], ncols = dm[2], crs = prj)
target[] <- 0
tf <- tempfile(fileext = ".tif")
raster::writeRaster(target, tf)
gt <- vapour::vapour_raster_info(tf)$geotransform

gt_dim_2ex <- function(x, dim) {
  xx <- c(x[1], x[1] + dim[1] * x[2]) - x[2]/2
  yy <- c(x[4] + dim[2] * x[6], x[4]) - x[6]/2
  raster::extent(c(xx, yy))
}

ex <- gt_dim_2ex(gt, dm)  ## extent from geotransform

f <- raadtools::topofile("etopo1")  ## local copy of a longlat raster of the globe
## warp the source grid to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(f, source_WKT = ll_wkt,target_geotransform = gt, target_dim = dm, target_WKT = wkt, band = 1)
#> setting projection
library(raster)
#> Loading required package: sp
plot(setValues(target, unlist(vals)))


system.time({
library(raster)
## 1. get a spec for the REST server image everyone wants
raw0 <- "https://raw.githubusercontent.com/OSGeo/gdal/master/gdal/frmts/wms/frmt_wms_arcgis_mapserver_tms.xml"
src <- c("World_Imagery", "World_Street_Map", "USA_Topo_Maps", "World_Terrain_Base",
         "Reference/World_Transportation",
         "Reference/World_Boundaries_and_Places_Alternate", "Reference/World_Boundaries_and_Places", "World_Shaded_Relief", 
         "World_Reference_Overlay", 
         "World_Physical_Map", "NatGeo_World_Map", "Ocean/World_Ocean_Base", 

         "Ocean/World_Ocean_Reference", "Canvas/World_Dark_Gray_Base", 
         "Canvas/World_Light_Gray_Base", "USA_Median_Household_Income", 
         "USA_Recent_Population_Change", "Elevation/World_Hillshade", 
         "Specialty/DeLorme_World_Base_Map", "Specialty/Soil_Survey_Map", "Specialty/World_Navigation_Charts")[1]
txt <- gsub("World_Street_Map", src, readLines(raw0))
#http://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer
filename <- tempfile(fileext = ".xml")
writeLines(txt, filename)

## 2. define an arbitrary projected area at a favourite region
## target coordinate system
## note I removed the false offsets
prj <-  "+proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
#prj <- "+proj=laea +lon_0=0 +lat_0=-45 +datum=WGS84"
#prj <- +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
## effectively a radius around my target projection centre
ex <- extent(-3e5, 4.5e5, -2.3e5, 7.8e5) + 1e5
wkt <- sf::st_crs(prj)$wkt

## 3. Create a target raster (we will populate its values from the warper)
par()
dm <- trunc(dev.size("px")*1)
target <- raster::raster(ex, nrows = dm[1], ncols = dm[2], crs = prj)
target[] <- 0

## 4. information about the source
## source coordinate system (here we don't need - GDAL often can't get one from a source, and needs user-input)
src_wkt <- ""

## 5. information the warper needs
#remotes::install_github("hypertidy/affinity")
gt <- affinity:::raster_to_gt(target)

## warp the source grid to the target grid
rgb_vals <- lapply(1:3, function(.x) {
  dirigible:::warp_in_memory_gdal_cpp(filename, source_WKT = src_wkt,
                                      target_geotransform = gt, 
                                      target_dim = dm[2:1], 
                                      target_WKT = wkt, 
                                      band = .x)
})

plotRGB(setValues(brick(target, target, target), do.call(cbind, unlist(rgb_vals, recursive = FALSE))), maxpixels = ncell(target))
})
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
#>  but +towgs84= values preserved

#>    user  system elapsed 
#>   1.999   0.036   3.512

^{Created on 2020-09-25 by the reprex package (v0.3.0)}

mdsumner commented 3 years ago

next iteration,

convenience for online source with gdalwebsrv
examples for building stars or raster
wkt and warp from dirigible
geotransform from affinity


## 1. get the terrain tms .xml
filename <- gdalwebsrv::lerc_file() # remotes::install_github("mdsumner/gdalwebsrv")

## 2. define an extent, in any projection
pt <- cbind(144.9631, -37.8136) + c(0, 0.1)
prj <-  "+proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
## get the projection as wkt (via GDAL)
wkt <- dirigible:::proj_to_wkt_gdal_cpp(prj)
xylim <- rep(reproj::reproj(pt, wkt, source = 4326)[,1:2], each = 2) + c(-1, 1, -1, 1) * 25000
(ex <- raster::extent(xylim))

## 3. Create a target raster (we will populate its values from the warper)
par() ## we use the device to get an appropriate pixel dims
dm <- trunc(dev.size("px")*1)
## spatial object, just a convenience for its extent,dim,crs (later we populate it for vis)
target <- raster::raster(ex, nrows = dm[1], ncols = dm[2], crs = prj)

## 4. information about the source
## source coordinate system (here we don't need - but we might if GDAL needs user-input)
src_wkt <- ""

## 5. information the warper needs
gt <- affinity:::raster_to_gt(target)  # remotes::install_github("hypertidy/affinity")
## warp the source grid to the target grid
vals <- dirigible:::warp_in_memory_gdal_cpp(filename, source_WKT = src_wkt,
                                            target_geotransform = gt, 
                                            target_dim = dm[2:1], 
                                            target_WKT = wkt, 
                                            band = 1)

## 6. set up and map the thing
## just convenience for value scaling to image
break_maker <- function(x, n = 10) {
  x <- unlist(x)
  quantile(x, seq(0, 1, length = n), na.rm = TRUE)
}

## create a raster from our template above
r <- raster::setValues(target, unlist(vals))

# library(stars)
# s <- stars::st_as_stars(st_bbox(target), nx = dm[2], ny = dm[1])
# s[[1]][] <- unlist(vals)
# image(s, col = gray.colors(25), breaks = break_maker(vals, n = 26));plot(lga, add = TRUE, col = NA, border = "black")

## make the map
library(ozmaps)
lga <- sf::st_crop(sf::st_transform(abs_lga, wkt), sf::st_bbox(target))
library(raster)
graphics.off()
par(xaxs = "i", yaxs = "i", mar = rep(0, 4))
image(r, col = gray.colors(25), breaks = break_maker(vals, n = 26), asp = 1, axes = FALSE)
plot(lga, add = TRUE, col = NA, border = "black")

dis-organization / dirigible

the in-memory raster warper #4