Reprojecting maps

Specific Domains Geo
1

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)

image
image802×642 168 KB

Thanks!

2

See GMT.jl

3

Thanks! Could you provide an example? From the docs I see mainly the plotting. Thanks in advance!

4

Rasters.reproject is the function in Rasters.jl, apologies I dont have time to write an example right now.

But have you tried that? The docs should specify the arguments you need.

5

Thanks, no had only tried Rasters.resample, but I could not get it to work.

Rasters.reproject gives this

julia> Rasters.reproject(ProjString("+proj=robin +datum=WGS84"), ras)
ERROR: StackOverflowError:
Stacktrace:
 [1] reproject(source::Nothing, target::ProjString, dim::X{Colon}, vals::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}) (repeats 79984 times)
   @ Rasters C:\Users\<username>\.julia\packages\Rasters\wUuLU\src\methods\reproject.jl:54

where ras was the return of world from the countryRaster function above.

Any idea?
Thanks,
Markus

6

Arrg sorry i was away from a computer resample is what you want:

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, missingval=0, fill=1, shape=:line)
crs = ProjString("+proj=robin +datum=WGS84")
wsg84_world = resample(world, res; crs)
Makie.heatmap(wsg84_world)

world
world960×1053 92.8 KB

I do want to make the syntax a little easier to understand for resample.

In future it will be better to use ArchGDAL.jl or Proj.jl to reproject the shapefile instead of resampling the raster, but my 2 PRs to make both of those things work on any geometries (like Shapefile.jl objects) are not quite finished.

1 Like
7

Thanks, yes in the real example I have a raster in the first place. Just used the county borders as a transparent MWE.

Unfortunately, I cannot reproduce, but get the error (Trace cut at [6]). Also I note that your output does not seem to be reprojected to Robinson.

julia> wsg84_world = Rasters.resample(world, reso; crs)
ERROR: BoundsError: attempt to access 0-element Vector{String} at index [1]
Stacktrace:
  [1] getindex
    @ .\array.jl:924 [inlined]
  [2] first
    @ .\abstractarray.jl:404 [inlined]
  [3] metadata(::ArchGDAL.RasterDataset{Int64, ArchGDAL.Dataset})
    @ Rasters C:\Users\mreichstein\.julia\packages\Rasters\wUuLU\src\sources\gdal.jl:231
  [4] dims(raster::ArchGDAL.RasterDataset{Int64, ArchGDAL.Dataset}, crs::Nothing, mappedcrs::Nothing)
    @ Rasters C:\Users\mreichstein\.julia\packages\Rasters\wUuLU\src\sources\gdal.jl:157
  [5] dims
    @ C:\Users\mreichstein\.julia\packages\Rasters\wUuLU\src\sources\gdal.jl:142 [inlined]
  [6] copyto!(dest::ArchGDAL.RasterDataset{Int64, ArchGDAL.Dataset}, bc::Base.Broadcast.Broadcasted{DimensionalData.DimensionalStyle{Base.Broadcast.DefaultArrayStyle{2}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(identity), Tuple{Raster{Int64, 2, Tuple{X{Projected{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, 
DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.NoMetadata, Nothing, Nothing, X{Colon}}}, Y{Projected{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.NoMetadata, Nothing, Nothing, Y{Colon}}}}, Tuple{}, Matrix{Int64}, Nothing, DimensionalData.Dimensions.LookupArrays.NoMetadata, Int64}}})
    @ DimensionalData C:\Users\mreichstein\.julia\packages\DimensionalData\1wSpb\src\array\broadcast.jl:48
8

Sorry I’m not at the computer. But maybe the example in the front page of the GMT.jl site

coast(region=:global, proj=:Winkel, frame=:g, area=10000,
      land=:burlywood4, water=:wheat1, show=true)

image

9

Yeah, it looks like there may be a bug in resample somewhere, but we may also be on different versions of some other package for you to get that error.

Can you make a Github issue? That’s usually the best thing to do when you get an error with Rasters.jl.

I’ll be away for a few days but having the issue there means I will get to it at some point after that.

10

FWIW, with some additional lines/tweaks, this is working now:

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(last, countries; to = Extents.Extent(X = (-180, 180), Y = (-90, 90)), res, missingval=0, fill=1, shape=:line)
world=setcrs(world, "+proj=longlat")
bounds = X(-180 .. 180), Y(-89.999 .. 89.999)
world2=crop(world; to=world[bounds...])
crs = ProjString("+proj=robin")
robin_world = resample(world2, 25000; crs, method=:near) # Define output resolution in new coord system (here 25000 m)
# Alternative:
#robin_world = resample(world2; crs, method=:near, size=(1440, 720)) # or define output size
heatmap(robin_world)
2 Likes