another python vs R example (WIP)

[mdsumner]

this was way more involved than I thought, I think my rasterio is too old

https://gist.github.com/mdsumner/8f070bc0bda18bb3722905f9afa1024b

work in progress ... but we should be able to get the same result as the R example, with about the same amount of code

[mdsumner]

here's a NASADEM example, taken from this notebook using OpenDataCube

https://github.com/opendatacube/cube-in-a-box/blob/main/notebooks/NASADEM.ipynb

(the code for the spatial region is out of sync with the actual data in the image)

do this in R with the following, it's not hillshade but I might explore that

library(vapour)
library(dsn) ## pak::pkg_install(c("hypertidy/dsn", "hypertidy/ximage"))
library(ximage)

## note that opendatacube expects res in top-down order 
##d <- gdal_raster_data(dsn::nasadem(), target_ex = c(146.0, 148.0, -43.5, -41.5), target_res = 0.001)
## I'd rather use the dimensions of the device than guess at a resolution (but this needs aspect ratio check too)
dm <- dev.size("px")
dm[which.max(dm)] <- 0
d <- gdal_raster_data(dsn::nasadem(), target_ex = c(146.0, 148.0, -43.5, -41.5), target_dim = dm)

## this just me finding how to get this particular colourmap from python
## https://stackoverflow.com/questions/18265935/how-do-i-create-a-list-with-numbers-between-two-values
# from matplotlib import cm  
# cm.gist_earth
# list2 = [x*(1/255) for x in range(0, 255)]
# cm.gist_earth(list2, bytes =True)

#cm <- reticulate::import("matplotlib.cm")
#colarr <- cm$gist_earth(seq(0, 1, length.out = 256), bytes = F)
#colnames(colarr) <- c("red", "green", "blue", "alpha")
#cols <- do.call(rgb, as.data.frame(colarr))
ximage::ximage(d, col = cols, asp = 1/cos(mean(attr(d, "extent")[3:4]) * pi/180))

Here's the cols saved so we don't need python call

cols <- c("#000000FF", "#01002BFF", "#010039FF", "#020043FF", "#03004EFF", 
"#030059FF", "#040064FF", "#05006EFF", "#050274FF", "#060574FF", 
"#070774FF", "#070975FF", "#080B75FF", "#090E75FF", "#091075FF", 
"#0A1275FF", "#0B1575FF", "#0B1776FF", "#0C1976FF", "#0D1C76FF", 
"#0D1E76FF", "#0E2076FF", "#0F2277FF", "#0F2577FF", "#102777FF", 
"#112977FF", "#112C77FF", "#122E77FF", "#133078FF", "#133278FF", 
"#143478FF", "#153778FF", "#153978FF", "#163B79FF", "#173D79FF", 
"#173F79FF", "#184179FF", "#194379FF", "#194579FF", "#1A477AFF", 
"#1B497AFF", "#1B4B7AFF", "#1C4D7AFF", "#1D4F7AFF", "#1D517BFF", 
"#1E537BFF", "#1F557BFF", "#1F577BFF", "#20597BFF", "#215A7BFF", 
"#215C7CFF", "#225E7CFF", "#23607CFF", "#23627CFF", "#24647CFF", 
"#25657DFF", "#25677DFF", "#26687DFF", "#276A7DFF", "#276B7DFF", 
"#286D7DFF", "#296F7EFF", "#29707EFF", "#2A727EFF", "#2B737EFF", 
"#2B757EFF", "#2C777FFF", "#2D787FFF", "#2D7A7FFF", "#2E7B7FFF", 
"#2F7D7FFF", "#2F7E7FFF", "#308080FF", "#30817EFF", "#31817DFF", 
"#31827CFF", "#32827BFF", "#328379FF", "#338378FF", "#338477FF", 
"#348576FF", "#348575FF", "#358673FF", "#358672FF", "#368771FF", 
"#368770FF", "#37886EFF", "#37896DFF", "#38896CFF", "#388A6BFF", 
"#388A6AFF", "#398B68FF", "#398C67FF", "#3A8C66FF", "#3A8D65FF", 
"#3B8D63FF", "#3B8E62FF", "#3C8E61FF", "#3C8F60FF", "#3D905FFF", 
"#3D905DFF", "#3E915CFF", "#3E915BFF", "#3F925AFF", "#3F9258FF", 
"#409357FF", "#409456FF", "#409455FF", "#419554FF", "#419552FF", 
"#429651FF", "#429650FF", "#43974FFF", "#43984DFF", "#44984CFF", 
"#44994BFF", "#45994AFF", "#459A49FF", "#479A47FF", "#4A9B46FF", 
"#4C9C47FF", "#4E9C47FF", "#509D48FF", "#539D48FF", "#559E49FF", 
"#579F4AFF", "#599F4AFF", "#5BA04BFF", "#5EA04BFF", "#60A14CFF", 
"#62A14CFF", "#64A24DFF", "#67A34EFF", "#69A34EFF", "#6BA44FFF", 
"#6DA44FFF", "#6FA450FF", "#72A551FF", "#74A551FF", "#76A652FF", 
"#78A652FF", "#7AA752FF", "#7CA753FF", "#7DA853FF", "#7FA853FF", 
"#81A854FF", "#82A954FF", "#84A954FF", "#86AA55FF", "#87AA55FF", 
"#89AB55FF", "#8AAB56FF", "#8CAB56FF", "#8EAC56FF", "#8FAC56FF", 
"#91AD57FF", "#93AD57FF", "#94AE57FF", "#96AE58FF", "#98AF58FF", 
"#99AF58FF", "#9BAF59FF", "#9DB059FF", "#9EB059FF", "#A0B15AFF", 
"#A2B15AFF", "#A3B25AFF", "#A5B25BFF", "#A7B25BFF", "#A8B35BFF", 
"#AAB35CFF", "#ABB45CFF", "#ADB45CFF", "#AFB55CFF", "#B0B55DFF", 
"#B2B65DFF", "#B4B65DFF", "#B5B65EFF", "#B7B75EFF", "#B7B65EFF", 
"#B8B55FFF", "#B8B45FFF", "#B9B35FFF", "#B9B360FF", "#B9B260FF", 
"#BAB160FF", "#BAB061FF", "#BBAF61FF", "#BBAE61FF", "#BCAD62FF", 
"#BCAC62FF", "#BCAC62FF", "#BDAB62FF", "#BDAA63FF", "#BEA963FF", 
"#BEA863FF", "#BFA764FF", "#BFA664FF", "#BFA664FF", "#C0A565FF", 
"#C0A465FF", "#C1A367FF", "#C2A46AFF", "#C3A46CFF", "#C4A46FFF", 
"#C5A571FF", "#C6A574FF", "#C7A676FF", "#C8A779FF", "#CAA87BFF", 
"#CBA97EFF", "#CCAA80FF", "#CDAA82FF", "#CEAB85FF", "#CFAC87FF", 
"#D0AD8AFF", "#D1AE8CFF", "#D3AF8FFF", "#D4B091FF", "#D5B194FF", 
"#D6B196FF", "#D7B298FF", "#D8B49BFF", "#D9B59DFF", "#DAB7A0FF", 
"#DBB8A2FF", "#DDB9A5FF", "#DEBBA7FF", "#DFBCAAFF", "#E0BDADFF", 
"#E1BFB0FF", "#E2C1B3FF", "#E3C3B6FF", "#E4C5B9FF", "#E6C7BCFF", 
"#E7C9BFFF", "#E8CBC2FF", "#E9CDC5FF", "#EACFC8FF", "#EBD1CBFF", 
"#ECD3CEFF", "#EDD5D1FF", "#EED8D4FF", "#F0DAD7FF", "#F1DCDAFF", 
"#F2DEDDFF", "#F3E1E0FF", "#F4E4E3FF", "#F5E7E6FF", "#F6EAE9FF", 
"#F7ECECFF", "#F9EFEFFF", "#FAF2F2FF", "#FBF5F5FF", "#FCF8F8FF", 
"#FDFBFBFF")

[mdsumner]

ok, so this was not too bad

library(vapour)
library(dsn)
library(ximage)

## note that opendatacube expects res in top-down order 
##d <- gdal_raster_data(dsn::nasadem(), target_ex = c(146.0, 148.0, -43.5, -41.5), target_res = 0.001)
## I'd rather use the dimensions of the device than guess at a resolution (but this needs aspect ratio check too)
dm <- dev.size("px")
dm[which.max(dm)] <- 0
d <- gdal_raster_data(dsn::nasadem(), target_ex = c(146.0, 148.0, -43.5, -41.5), target_dim = dm)

## this just me finding how to get this particular colourmap from python
## https://stackoverflow.com/questions/18265935/how-do-i-create-a-list-with-numbers-between-two-values
# from matplotlib import cm  
# cm.gist_earth
# list2 = [x*(1/255) for x in range(0, 255)]
# cm.gist_earth(list2, bytes =True)

## the data
m <- matrix(d[[1L]], attr(d, "dimension")[2], byrow = TRUE)
rs <- rayshader::lamb_shade(m, zscale = 0.005 )

cm <- reticulate::import("matplotlib.cm")
colarr <- cm$gist_earth(seq(0, 1, length.out = 256), bytes = F)
colnames(colarr) <- c("red", "green", "blue", "alpha")
cols <- do.call(rgb, as.data.frame(colarr))
par(mfrow = c(1, 2))

asp <- 1/cos(mean(attr(d, "extent")[3:4]) * pi/180)
ximage::ximage(rs[nrow(rs):1, ], col = grey.colors(256), extent = attr(d, "extent"), 
               asp = asp)

ximage::ximage(d, attr(d, "extent"), col = scales::alpha(cols, 0.85), add = TRUE)
ximage::ximage(d, attr(d, "extent"), col =cols, add = FALSE, asp = asp)

hillshaded version on the left

[mdsumner]

I got some examples to work with sf, afaik you can't use these cogs or webservers reliably with stars at all. The only way that works is to redirect to a file with "sf::gdal_utils('warp', ...".

the sf gdal_utils warper doesn't work well with these huge sources (it would be ok for a very small region).

the sf gdal_utils warp can work if you set the options as per gdalwarp (you need at least -tr, -ts, or -te and for te it should be a small region unless also using ts or tr)

## pak::pkg_install("hypertidy/dsn")
src_elev <- dsn::gebco()
src_rgb <- dsn::wms_arcgis_mapserver_ESRI.WorldImagery_tms()

library(sf)

sf::gdal_utils("warp", src_elev, tf <- tempfile(fileext = ".tif"), 
  options = c("-ts", 1024, 512))
graphics.off(); plot(stars::read_stars(tf))

sf::gdal_utils("warp", src_rgb, tf <- tempfile(fileext = ".tif"), 
  options = c("-ts", 1024, 1024, "-t_srs", "+proj=laea"))
graphics.off(); plot(stars::read_stars(tf), rgb = 1:3)

[joshualerickson]

I had no idea sf had a gdal abstraction!? That function looks sweet. Need to start adding these to the repo…

[mdsumner]

sf has all the gdal stuff for r-spatial/ (stars uses it, for example)

terra has all the gdal stuff for rspatial/ :)