MakieOrg / GeoMakie.jl

Geographical plotting utilities for Makie.jl
https://geo.makie.org
MIT License
173 stars 25 forks source link

Problems with cropping subplot lat/lon extents #122

Open JordiBolibar opened 2 years ago

JordiBolibar commented 2 years ago

I'm new to GeoMakie, and I'm trying to start using it as my main tool for plotting maps in Julia. I've been struggling to implement a plot with multiple subplots with different lat/lon ranges. Here's the code (not reproducible):

theme = Theme(fontsize=18, font="TeX Gyre Heros")
set_theme!(theme)

# Ice surface velocities
figV = Makie.Figure(resolution=(1200, 1800))
n = 5
m = floor(Int, length(gdirs_climate[2])/n)

axsV = []
k = 1
for i in 1:n
    for j in 1:m
        push!(axsV, GeoAxis(figV[i, j], source = "+proj=tmerc +lon_0=$(gdirs[k].cenlon) +lat_0=$(gdirs[k].cenlat)", 
                    dest = "+proj=tmerc +lon_0=$(gdirs[k].cenlon) +lat_0=$(gdirs[k].cenlat)", 
                    lonlims=(gdirs[k].extent_ll[1,1], gdirs[k].extent_ll[1,2]), 
                    latlims=(gdirs[k].extent_ll[2,1], gdirs[k].extent_ll[2,2]), 
                    coastlines = true, remove_overlapping_ticks = true))
        k += 1
    end
end

tightlimits!.(axsV)
Vs, hm_Vs = [],[]
for (i, ax) in enumerate(axsV)
    ax.aspect = DataAspect()
    rgi_id = gdirs[i].rgi_id
    name = gdirs[i].name
    ax.title = "$name - $rgi_id \n $(gdirs[i].cenlat)° - $(gdirs[i].cenlon)°"
    Vx = PDE_refs["V̄x_refs"][i]
    Vy = PDE_refs["V̄y_refs"][i]
    V = (abs.(Vx) .+ abs.(Vy))./2
    lons = gdirs[i].extent_ll[1,:]
    lats = gdirs[i].extent_ll[2,:]

    push!(Vs,V)
    lonrange = collect(range(lons[1], lons[2], length=size(V)[1]))
    latrange = collect(range(lats[1], lats[2], length=size(V)[2]))
    push!(hm_Vs, Makie.heatmap!(ax, lonrange, latrange, fillZeros(V), colormap=:speed))
end
minV = minimum(minimum.(Vs))
maxV = maximum(maximum.(Vs)) 
Colorbar(figV[2:3,m+1], limits=(minV,maxV), colormap=:speed, label="Ice surface velocity")
Label(figV[0, :], text = "Glacier dataset", textsize = 30)
if display
    display(figV)
end
Makie.save(joinpath(root_plots, "glaciers_V.pdf"), figV, pt_per_unit = 1)

Which produces the following plot: image

What am I doing wrong? How can I correctly crop the subplot extent in order to "zoom" on the given lat/lon extent? Sorry for the noob question, but I haven't been able to figure this out with the documentation. Thanks in advance!

asinghvi17 commented 2 years ago

I suspect that tightlimits! might be doing something wrong...maybe try datalims! instead, which is GeoMakie's version of the function?

Otherwise I'm kind of tied up now but can test this out in a bit.

asinghvi17 commented 5 months ago

I still occasionally see these kinds of errors with interactive rectangle zoom, but never with zooming by setting limits.