Plotting data with 2D lat/lons in geomakie

k945 · October 13, 2022, 8:09am

Is it possible to plot data with 2d lat/lons using GeoMakie? I have so far been unable to properly project WRF simulation data onto a map.

Using test data, this is what I have tried so far:

using GeoMakie, CairoMakie, NetCDF

mu = "test.nc"

# source projection gotten earlier from wrfout file
source = "+proj=ob_tran +o_proj=latlon +o_lon_p=-141.38 +o_lat_p=-60.31 +lon_0=147.63 +R=6370000 +nadgrids=@null +no_defs +type=crs"

# Dest = cartopy.crs.PlateCarree().to_proj4()  (as has been uused in https://wrf-python.readthedocs.io/en/latest/plot.html).
dest = "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +to_meter=111319.490793274 +no_defs +type=crs"


landmask = ncread(mu, "landmask")
lat = ncread(mu, "lat")
lon = ncread(mu, "lon")
lon = map(x->x<0 ? 360 + x : x, lon) # It seems that 

begin 
    fig = Figure()

#  better projection when setting dest = source rather than dest = dest
    ax = GeoAxis(fig[1,1]; coastlines = true, source = source, dest = source) 
    sf = heatmap!(ax, lon,lat, landmask, colormap = :Hiroshige, colorrange = (0,1))  
    datalims!(ax) 
    fig
end

If it is possible, I would rather do this in GeoMakie instead of using PyPlot (like I had previously done).

Raf · October 13, 2022, 8:18am

You may be missing that NetCDF can have a source projection different to the values in the dimension variables - because they are a vector rather than a step-range or geotransform. You may need to convert those values from epsg:4326 to your source projection first (or something).

Rasters.jl can handle these conversions automatically, but isn’t integrated with GeoMakie.jl yet. We should fix that, you should just be able to plot any raster like this without doing this work.

(You can also just do this plot with Rasters.jl + Plots.jl and it should be fine plot(Raster(filename)))

k945 · October 13, 2022, 5:43pm

I got the source projection using the wrf-python get_cartopy function. It does not look like any conversion/transformation is taking place in the following plotting example, so I am assuming that the source projection and dimension variables correspond to each other.

Unfortunately, I am have not been able to use Rasters.jl because it does not seem to be able to load WRF (netcdf) data that consists of 2D lat/lons. For example,

julia> Rasters("test.nc", key = :landmask)

julia> Raster(mu; key= :landmask)
typeof(index) = Matrix{Float32}
ERROR: StackOverflowError:
Stacktrace:
  [1] isrev(x::Type) (repeats 17293 times)
    @ DimensionalData.Dimensions.LookupArrays /data2/user/.julia/packages/DimensionalData/DkTAI/src/LookupArrays/lookup_traits.jl:70
  [2] _ncdspan(index::Matrix{Float32}, order::DimensionalData.Dimensions.LookupArrays.Unordered)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:327
  [3] _ncdlookup(ds::NCDataset{Nothing}, dimname::String, D::Type, index::Matrix{Float32}, metadata::DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, crs::EPSG, mappedcrs::EPSG)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:290
  [4] _ncdlookup(ds::NCDataset{Nothing}, dimname::String, D::Type, crs::EPSG, mappedcrs::EPSG)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:264
  [5] _ncddim(ds::NCDataset{Nothing}, dimname::String, crs::EPSG, mappedcrs::EPSG)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:231
  [6] #725
    @ /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:147 [inlined]
  [7] map
    @ ./tuple.jl:222 [inlined]
  [8] dims
    @ /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:146 [inlined]
  [9] Raster(ds::NCDatasets.CFVariable{Float32, 2, NCDatasets.Variable{Float32, 2, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Nothing, Tuple{}, Nothing, Nothing, Nothing, Nothing, Nothing}}}, filename::String, key::Symbol; crs::Nothing, mappedcrs::Nothing, dims::Nothing, refdims::Tuple{}, name::Symbol, metadata::DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.NCDfile, Dict{Symbol, Any}}, missingval::Missing, source::Type, write::Bool, lazy::Bool)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/array.jl:273
 [10] Raster(ds::NCDatasets.CFVariable{Float32, 2, NCDatasets.Variable{Float32, 2, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Nothing, Tuple{}, Nothing, Nothing, Nothing, Nothing, Nothing}}}, filename::String, key::Symbol)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/array.jl:264
 [11] Raster(ds::NCDataset{Nothing}, filename::String, key::Symbol; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:34
 [12] Raster(ds::NCDataset{Nothing}, filename::String, key::Symbol)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:32
 [13] (::Rasters.var"#51#52"{Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, String})(ds::NCDataset{Nothing})
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/array.jl:261
 [14] _open(f::Rasters.var"#51#52"{Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, String}, ::Type{Rasters.NCDfile}, ds::NCDataset{Nothing}; key::Nothing, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:222
 [15] _open
    @ /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:221 [inlined]
 [16] #745
    @ /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:218 [inlined]
 [17] NCDataset(::Rasters.var"#745#746"{Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Rasters.var"#51#52"{Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, String}}, ::String, ::Vararg{String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NCDatasets /data2/user/.julia/packages/NCDatasets/ipGBH/src/dataset.jl:241
 [18] NCDataset(::Function, ::String, ::Vararg{String})
    @ NCDatasets /data2/user/.julia/packages/NCDatasets/ipGBH/src/dataset.jl:238
 [19] #_open#744
    @ /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:217 [inlined]
 [20] _open(f::Function, ::Type{Rasters.NCDfile}, filename::String)
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/sources/ncdatasets.jl:214
 [21] _open(f::Function, filename::String; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Rasters /data2/user/.julia/packages/Rasters/fANS1/src/convenience.jl:20
 [22] _open
    @ /data2/user/.julia/packages/Rasters/fANS1/src/convenience.jl:19 [inlined]
 [23] #Raster#50
    @ /data2/user/.julia/packages/Rasters/fANS1/src/array.jl:259 [inlined]
 [24] top-level scope
    @ REPL[39]:1

For now, I will revert back to using my old workflow in order to plot the data. Thanks for the help.

On an unrelated matter: Thanks to you and everyone else that has worked on DimensionalData.jl. It has been very useful.

Raf · October 13, 2022, 6:47pm

Interesting - yes I haven’t implemented handling of matrix dimension variables yet, not being sure how common it is.

But please make github issues whenever you hit errors like this. I really mean it, they are the main gauge I have of things people need and what to prioritise - basically if there is no issue, there will be no feature unless I need it myself. There are too many things to do to implement things that people don’t want and I don’t need.

(With the plot I still think the issue will be the crs - those matrix variables are usually a translation layer from the underlying projection to another projection - often epsg:4326)

k945 · November 2, 2022, 5:08am

Which repo would you prefer I open up the issue in: DimensionalData or Rasters?

k945 · November 2, 2022, 5:57am

It appears that the 2d lat/lon data is valid wgs84.

I am not sure why this is the case but plotting using surface! instead of heatmap! seems to plot line up with the with the map.

begin 
    source = "+proj=eqc +o_alpha=0 +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +to_meter=111319.490793274 +no_defs +type=crs"
    dest = source

    fig = GeoMakie.Figure()
    ax = GeoAxis(fig[1,1]; source = source, dest = dest, coastlines = true,  lonlims = automatic, latlims = automatic)
 
    sf = surface!(ax, lon, lat, landmask[:,:], colormap = :viridis, shading = false)
    Colorbar(fig[2,1], sf; vertical = false)

    datalims!(ax) 

    fig
end

The question now is, what is an appropriate transformation? For example, setting the destination

 dest = "+proj=ob_tran +o_proj=latlon +o_lon_p=-141.380004882812 +o_lat_p=-60.310001 +lon_0=147.630001068115 +R=6370000 +nadgrids=@null +no_defs +type=crs"

produces a flipped map

Which is not quite what I want.

lazarusA · November 2, 2022, 6:00am

quick note: Using surface is preferable to heatmap, which is does not work in most cases, but surface does. GeoMakie.jl · GeoMakie.jl

k945 · November 2, 2022, 6:28am

The documentation (at least for me) is not very clear with regards to this. It states that

With CairoMakie, we recommend that you use image!(ga, ...) or heatmap!(ga, ...) to plot images or scalar fields into ga::GeoAxis.

It then goes on to say

with GLMakie, …, these methods do not work

The way you have phrased it here is much clearer.

I will keep this in mind. Thanks

Raf · November 2, 2022, 9:37am

Generally either is fine, I can move them.

This probably needs some work in both, so Rasters.jl.