Dear community,
does anyone have a recipe/example of how to reproject a map, e.g. from lat/lon to Robinson (for my use case just global is enough)? Please note I do not just want to plot (cf. GeoMakie), but get a Raster or matrix back.
Essentially in the MWE below I would like to avoid using RCall:
using GLMakie
using RCall
using Rasters
import Downloads, Shapefile
@rlibrary terra
function countryRaster(; res=0.5)
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "country_borders.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)
countries = Shapefile.Handle(shapefile_name).shapes
world = Rasters.rasterize(countries; to = Extents.Extent(X = (-180, 180), Y = (-90, 90)), res=res, missingval=0, fill=1, shape=:line)
return world.data.data
end
function proj(grid, outcrs="+proj=robin +datum=WGS84")# +units=m +no_defs")
#grid = permutedims(grid, (2,1))
jras = R"""
terra::rast($(permutedims(grid, (2,1))), crs = "+proj=longlat +datum=WGS84 +no_defs", extent=ext(rast())) %>% terra::project($outcrs) %>% as.matrix(wide=T)
""" |> rcopy
return permutedims(jras, (2,1))
end
heatmap(countryRaster() |> proj)
Thanks!
