# Plotting data with 2D lat/lons in geomakie

**URL:** <https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648>\
**Category:** Geo\
**Tags:** question, package, plotting\
**Created:** [October 13, 2022, 8:09am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648 "2022-10-13T08:09:08Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![k945](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/k945/32/8906_2.png) [@k945](https://discourse.julialang.org/u/k945)\
**Post date:** [October 13, 2022, 8:09am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/1 "2022-10-13T08:09:09Z")

</div>

Is it possible to plot data with 2d lat/lons using GeoMakie? I have so far been unable to properly project [WRF](https://www.mmm.ucar.edu/models/wrf) simulation data onto a map.

 ![test](https://global.discourse-cdn.com/julialang/original/3X/2/6/26758e3c277ab60806fd992d1d0130b2c7966cde.png)  
Using [test data](https://cloudstor.aarnet.edu.au/plus/s/QuzKe9tg8YIfcZF), this is what I have tried so far:

```julia
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](https://discourse.julialang.org/t/pyplot-not-projecting-netcdf-as-expected/9446)).

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 13, 2022, 8:18am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/2 "2022-10-13T08:18:22Z")

</div>

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))`)

---

<div class="post-metadata">

**Author:** ![k945](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/k945/32/8906_2.png) [@k945](https://discourse.julialang.org/u/k945)\
**Post date:** [October 13, 2022, 5:43pm UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/3 "2022-10-13T17:43:11Z")

</div>

I got the source projection using the [wrf-python](https://wrf-python.readthedocs.io) get\_cartopy function. It does not look like any conversion/transformation is taking place in the following [plotting example](https://wrf-python.readthedocs.io/en/latest/plot.html#matplotlib-with-cartopy), 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
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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 13, 2022, 6:47pm UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/4 "2022-10-13T18:47:48Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![k945](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/k945/32/8906_2.png) [@k945](https://discourse.julialang.org/u/k945)\
**Post date:** [November 2, 2022, 5:08am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/5 "2022-11-02T05:08:39Z")

</div>

> [@Raf](#):
>
> But please make github issues whenever you hit errors like this

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

---

<div class="post-metadata">

**Author:** ![k945](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/k945/32/8906_2.png) [@k945](https://discourse.julialang.org/u/k945)\
**Post date:** [November 2, 2022, 5:57am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/6 "2022-11-02T05:57:50Z")

</div>

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.

```julia
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

```

 ![plot](https://global.discourse-cdn.com/julialang/original/3X/4/f/4f90e8be100f4af2e1a38415989404597c02eb80.png)

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

```julia
 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

 ![malformed](https://global.discourse-cdn.com/julialang/original/3X/e/c/ec62d1ac7095d920267f2d293068e6c09a4922b6.png)

Which is not quite what I want.

---

<div class="post-metadata">

**Author:** ![lazarusA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lazarusa/32/6571_2.png) [@lazarusA](https://discourse.julialang.org/u/lazarusA)\
**Post date:** [November 2, 2022, 6:00am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/7 "2022-11-02T06:00:44Z")

</div>

quick note: Using surface is preferable to heatmap, which is does not work in most cases, but surface does. [GeoMakie.jl · GeoMakie.jl](https://geo.makie.org/stable/#Gotchas-2)

---

<div class="post-metadata">

**Author:** ![k945](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/k945/32/8906_2.png) [@k945](https://discourse.julialang.org/u/k945)\
**Post date:** [November 2, 2022, 6:28am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/8 "2022-11-02T06:28:23Z")

</div>

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.

> [@lazarusA](#):
>
> Using surface is preferable to heatmap, which is does not work in most cases, but surface does

I will keep this in mind. Thanks

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [November 2, 2022, 9:37am UTC](https://discourse.julialang.org/t/plotting-data-with-2d-lat-lons-in-geomakie/88648/9 "2022-11-02T09:37:24Z")

</div>

Generally either is fine, I can move them.

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