Skip to content

Instantly share code, notes, and snippets.

Avatar

Michael Sumner mdsumner

View GitHub Profile
View graticule_masking.md
  library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1
edge0 <- function(x, y, ndiscr = 18) {
  dx <- if(length(x) > 1) diff(x) else diff(y)
  sf::st_set_crs(sf::st_segmentize(sf::st_sfc(sf::st_linestring(cbind(x, y))), dx/ndiscr), 
               "OGC:CRS84")
}
north <- function(x = c(-180, 180), y = 90, ndiscr = 18) {
 edge0(x, y, ndiscr)  
View Adams_world_virt_earth.md
  library(sf)
#> Linking to GEOS 3.9.2, GDAL 3.3.2, PROJ 8.2.0
library(ggplot2)
library(rnaturalearth)

world_g <- st_graticule(ndiscr = 10000,
                        margin = 0.00000000001) |> 
  st_transform("+proj=adams_ws2") |> 
  st_geometry()
@mdsumner
mdsumner / wms_data.md
Last active Nov 12, 2021
Extracting data from WMS using R
View wms_data.md
View ice_50k.md
  library(raadtools)
rgdal::set_thin_PROJ6_warnings(TRUE)

files <- raadfiles::amsr2_daily_files()

library(dplyr)
files <- filter(files, between(as.Date(date), as.Date("2014-01-01"), as.Date("2021-01-01")))

## use as a template
View gdal_rest.md
gdal_translate "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json" wms.xml -of WMS
gdalwarp wms.xml -te -1e6 -1e6 1e6 1e6 -ts 512 512 -t_srs "+proj=laea +lon_0=147 +lat_0=-42" out.tif -of COG
library(terra)
plotRGB(terra::rast("out.tif"))
View terra_project_issue.md

Why does terra take forever (or never complete ...)?

## remotes::install_github("hypertidy/gdalio")
library(gdalio) 

## warp target 
extent <- c(-2e5, 2.4e5,  -2e5, 2.84e5)
dimension <- c(1024, 1024)
projection <- "+proj=lcc +lon_0=146 +lat_0=-42 +lat_1=-40 +lat_0=-44 +datum=WGS84"
View ggmap_360.md

get a range of ggmap that crosses longitude 180

library(ggmap)
mapwest <- get_map(location = c(left = -180, bottom = -40,right = -120,top = 20), source = "google", maptype = "terrain")
mapeast <- get_map(location = c(left = 110, bottom = -40,right = 180,top = 20), source = "google", maptype = "terrain")

## combine the images together (and then we'll reconstruct their class/metadata)
map <- as.raster(cbind(as.matrix(mapeast), as.matrix(mapwest)))
@mdsumner
mdsumner / ssh_reproject.md
Last active Oct 6, 2021
reproject SSH netcdf data to polar stereographic with {gdalio}
View ssh_reproject.md

Should be about 30x faster than using raster::projectRaster

    library(raster)

    library(gdalio)
    library(vapour)
    library(raadtools)
    
View mapbox_wmts.md

https://twitter.com/mdsumner/status/1440933887354441731?s=20

## build mapbox WMTS data source name from a custom style
mapbox_wmts <- function(type = "", key = NULL) {
  if (nchar(type) < 1) stop("type must be a useable mapbox style name e.g. 'styles/v1/<username>/<style_id>'")
  if (is.null(key)) {
    key <- Sys.getenv("MAPBOX_API_KEY")
    if (is.null(key)) stop("provide mapbox 'key' or via MAPBOX_API_KEY env var")
View opentopo_alos_warp.R
alos <- "/vsicurl/https://opentopography.s3.sdsc.edu/raster/AW3D30/AW3D30_global.vrt"
##eng_wales in Hugh's example (but I reduce the dimensions, for now)
g <- list(extent = c(82667.22, 655646.4, 5337.574, 657533.7), dimension = c(28575, 29006), projection = "EPSG:27700")
library(gdalio)
gdalio_set_default_grid(g)
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
r <- gdalio_raster(alos)