# GeoMakie: how to plot shapefile lines in the same GeoAxis as a heatmap

**URL:** <https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380>\
**Category:** New to Julia\
**Tags:** geo, glmakie, shapefile, netcdf, geomakie\
**Created:** [September 13, 2024, 5:35pm UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380 "2024-09-13T17:35:29Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![oscarvdvelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscarvdvelde/32/202157_2.png) [@oscarvdvelde](https://discourse.julialang.org/u/oscarvdvelde)\
**Post date:** [September 13, 2024, 5:35pm UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/1 "2024-09-13T17:35:30Z")

</div>

I am trying GeoMakie for a simple visualization of KNMI radar data (Netherlands), which comes in netCDF format. The netCDF provides grid values, the x and y axes and the proj4. I am able to successfully plot that georeferenced with:

```julia
ga = GeoAxis(fig[1,1], dest=radar[1].prj, source=radar[1].prj)
heatmap!(ga,radar[1].x, radar[1].y, radar[1].val)

```

where radar[1].prj is  
`"+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0 +units=km"`  
It plots a reference grid of latitude and longitude (despite not having a frame around it, which I would prefer).  
However, I am struggling to figure out the right way to plot a shapefile of country borders or the built in `GeoMakie.coastlines` correctly on top of it. Is there a good example available for this?

---

<div class="post-metadata">

**Author:** ![AMB](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amb/32/211738_2.png) [@AMB](https://discourse.julialang.org/u/AMB)\
**Post date:** [September 13, 2024, 6:43pm UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/2 "2024-09-13T18:43:35Z")

</div>

When it comes to drawing the coastlines, what have you tried?

```julia
lines!(ga, GeoMakie.coastlines()) 

```

---

<div class="post-metadata">

**Author:** ![oscarvdvelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscarvdvelde/32/202157_2.png) [@oscarvdvelde](https://discourse.julialang.org/u/oscarvdvelde)\
**Post date:** [September 13, 2024, 8:58pm UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/3 "2024-09-13T20:58:14Z")

</div>

That is indeed what I tried:

```julia
julia> ga = GeoAxis(fig[1,1], dest=radar[1].prj, source=radar[1].prj)
GeoAxis()

julia> heatmap!(ga,radar[1].x,radar[1].y,radar[1].val)
Heatmap{Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float32}}}
# correctly displays geo-referenced data

julia> countries = Shapefile.Handle(raw"C:\Users\oscar\Documents\Work\Julia\LMA\map\world-administrative-boundaries.shp").shapes
256-element Vector{Union{Missing, Shapefile.Polygon}}:
 Polygon(24 Points) # etc...

julia> poly!(ga,countries)
# or:
julia> lines!(ga,GeoMakie.coastlines())
# displays them small, all the way at the north pole

```

I have tried to play with the source and destination projection strings within the GeoAxis, before/after plotting.  
It looks like I cannot set `ga.source[] = "+proj=longlat"` after plotting the data and before adding the countries, because it also changes how the radar data are interpreted.  
Perhaps I should first convert the data grid to latitude longitude (how…). But when I do that, I also wouldn’t need GeoMakie anymore, as projections are not so critical for displaying a limited area up to a few hundred kilometers wide.

I have also tried to load the netCDF with `Rasters.jl`, as it may have functions for conversion, but it did not pick up the information properly:

```julia
julia> ras = Raster(radfile)
╭────────────────────────────────────────────────────────────╮
│ 0-dimensional Raster{Union{Missing, Float64},0} projection 
# only shows projection data...

```

---

<div class="post-metadata">

**Author:** ![AMB](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amb/32/211738_2.png) [@AMB](https://discourse.julialang.org/u/AMB)\
**Post date:** [September 14, 2024, 12:13am UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/4 "2024-09-14T00:13:33Z")

</div>

Maybe have a look at [GitHub - JuliaGeo/Proj.jl: Julia wrapper for the PROJ cartographic projections library](https://github.com/JuliaGeo/Proj.jl)

---

<div class="post-metadata">

**Author:** ![oscarvdvelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscarvdvelde/32/202157_2.png) [@oscarvdvelde](https://discourse.julialang.org/u/oscarvdvelde)\
**Post date:** [September 14, 2024, 2:46am UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/5 "2024-09-14T02:46:58Z")

</div>

Thanks, it is easily converted. This works now. Perhaps it could be worked into an example for the GeoMakie site.

```julia
trans = Proj.Transformation(prj,"WGS84",always_xy=true)
xy = [[x[i],y[j]] for i in 1:dims[1], j in 1:dims[2]]
lonlats = trans.(xy)
lon = reshape(getindex.(lonlats,1),dims)
lat = reshape(getindex.(lonlats,2),dims)

fig = nothing; fig = Figure(size = (900,900))
ga = GeoAxis(fig[1,1], limits=((0,10),(49,55)), dest="+proj=merc +lat_ts=52.0")
surface!(ga, lon,lat, values, colorrange=(0,percentile(skipmissing(values),99)))
poly!(ga, countries; strokecolor=(:gray,0.5), strokewidth = 1.25, overdraw = true)

```

The opposite transformation, from `countries` or `GeoMakie.coastlines()` WGS84 to the radar projection seems to need some additional conversion step before:

```julia
julia> trans.(GeoMakie.coastlines())
ERROR: ArgumentError: Argument is not a Point geometry

```

Could anyone help with that?

Update: corrected the alignment of x and y for the coordinate transform. Update 2: Succeeded plotting, here is the result!

 ![Screenshot 2024-09-14 073723](https://global.discourse-cdn.com/julialang/original/3X/3/e/3e83cb10428015539e1dc062e02f467885917ff0.jpeg)

---

<div class="post-metadata">

**Author:** ![AMB](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amb/32/211738_2.png) [@AMB](https://discourse.julialang.org/u/AMB)\
**Post date:** [September 15, 2024, 7:14pm UTC](https://discourse.julialang.org/t/geomakie-how-to-plot-shapefile-lines-in-the-same-geoaxis-as-a-heatmap/119380/6 "2024-09-15T19:14:03Z")

</div>

That looks amazing! 👍