Try using inverse transforms and basing limits off transformed space

src/geoaxis.jl

@@ -7,8 +7,9 @@ Makie.@Block GeoAxis begin
     targetlimits::Observable{Rect2d}
     # "Final limits in input space"
     finallimits::Observable{Rect2d}
-    # Final limits in transformed space
-    transformedlimits::Observable{Rect2d}
+    inputlimits::Observable{Union{Nothing, Rect2d}}
+    # The default transformation, to cache and save on calls to Proj!
+    transform_func::Observable{Any}
     # interaction stuff
     mouseeventhandle::Makie.MouseEventHandle
     scrollevents::Observable{Makie.ScrollEvent}
@@ -289,25 +290,23 @@ Makie.can_be_current_axis(::GeoAxis) = true
 
 function Makie.initialize_block!(axis::GeoAxis)
 
+    ptrans = create_transform(axis.source_projection, axis.target_projection)
+    setfield!(axis, :transform_func, ptrans)
+
     scene = axis_setup!(axis)
     setfield!(axis, :elements, Dict{Symbol,Any}())
 
-
-    ptrans = create_transform(axis.source_projection, axis.target_projection)
-
-    draw_geoaxis!(axis, axis.target_projection, axis.elements, false)
+    draw_geoaxis!(axis, axis.transform_func, axis.elements, false)
 
     return axis
 end
 
 # do the axis drawing
 
-function draw_geoaxis!(ax::GeoAxis, target_projection, elements, remove_overlapping_ticks)
+function draw_geoaxis!(ax::GeoAxis, transformation, elements, remove_overlapping_ticks)
     topscene = ax.blockscene
     scene = ax.scene
 
-    transformation = GeoMakie.create_transform(Observable("+proj=longlat +datum=WGS84"), target_projection)
-
     xgridpoints = Observable(Point2f[])
     ygridpoints = Observable(Point2f[])
 
@@ -333,7 +332,7 @@ function draw_geoaxis!(ax::GeoAxis, target_projection, elements, remove_overlapp
 
     # First we establish the spine points
 
-    lift(ax.finallimits, ax.xticks, ax.xtickformat, ax.yticks, ax.ytickformat, ax.xminorticks, ax.yminorticks, ax.scene.px_area, transformation, ax.npoints) do limits, xticks, xtickformat, yticks, ytickformat, xminor, yminor, pxarea, transform_func, npoints
+    lift(ax.inputlimits, ax.xticks, ax.xtickformat, ax.yticks, ax.ytickformat, ax.xminorticks, ax.yminorticks, ax.scene.px_area, transformation, ax.npoints) do limits, xticks, xtickformat, yticks, ytickformat, xminor, yminor, pxarea, transform_func, npoints
 
         lmin = minimum(limits)
         lmax = maximum(limits)
@@ -583,7 +582,7 @@ function draw_geoaxis!(ax::GeoAxis, target_projection, elements, remove_overlapp
     setproperty!.(values(elements), Ref(:yautolimits), Ref(false))
 
     # finally, make sure that lift runs again - for some reason, it doesn't work directly
-    notify(ax.transformedlimits)
+    notify(ax.inputlimits)
     notify(ax.finallimits)
 
     return nothing

src/makie-axis.jl

@@ -9,52 +9,69 @@ Base.convert(::Type{Rect2d}, x::Rect2) = Rect2d(x)
 # Makie.xautolimits(ga::GeoAxis) = (-180, 180)
 # Makie.yautolimits(ga::GeoAxis) = (-90, 90)
 
-function geodefaultlimits(::Tuple{Nothing, Nothing}, source_projection, target_projection)
-    return Rect2d(-180, -90, 360, 180)
+
+geodefaultlimits(axis::GeoAxis, limits) = geodefaultlimits(limits, axis.source_projection[], axis.target_projection[])
+
+# override for the automatic case
+function geodefaultlimits(axis::GeoAxis, ::Tuple{Nothing, Nothing})
+    # if isassigned(axis.scene)
+    #     transformed_data_limits = data_limits(axis.scene)
+    #     need_auto_limits = any(isinf.(origin(transformed_data_limits))) || any(isinf.(widths(transformed_data_limits)))
+    # else
+        need_auto_limits = true
+    # end
+
+    if need_auto_limits
+         # if the transformed data limits are infinite, nothing has been plotted
+        return Makie.apply_transform(axis.transform_func[], Rect2d(-180, -90, 360, 180)) # TODO: make this figure out what the appropriate limits are automagically
+    else
+        return transformed_data_limits
+    end
 end
 
-function geodefaultlimits(limits::Tuple{NTuple{2, <: Real}, Nothing}, source_projection, target_projection)
-    return Rect2d(limits[1][1], -90, limits[1][2] - limits[1][1], 180)
+function geodefaultlimits(axis::GeoAxis, limits::Tuple{NTuple{2, <: Real}, Nothing})
+    return Makie.apply_transform(axis.transform_func[], Rect2d(limits[1][1], -90, limits[1][2] - limits[1][1], 180))
 end
 
-function geodefaultlimits(limits::Tuple{ Nothing, NTuple{2, <: Real}}, source_projection, target_projection)
-    return Rect2d(-180, limits[2][1], 360, limits[2][2] - limits[2][1])
+function geodefaultlimits(axis::GeoAxis, limits::Tuple{ Nothing, NTuple{2, <: Real}})
+    return Makie.apply_transform(axis.transform_func[], Rect2d(-180, limits[2][1], 360, limits[2][2] - limits[2][1]))
 end
 
-function geodefaultlimits(limits::Tuple, source_projection, target_projection)
+function geodefaultlimits(axis::GeoAxis, limits::Tuple)
     (xmin, xmax), (ymin, ymax) = limits
-    return Rect2d(xmin, ymin, xmax - xmin, ymax - ymin)
+    return Makie.apply_transform(axis.transform_func[], Rect2d(xmin, ymin, xmax - xmin, ymax - ymin))
 end
 
-function geodefaultlimits(limits::Rect{2, <: Real}, source_projection, target_projection)
-    return Rect2d(limits)
+function geodefaultlimits(axis::GeoAxis, limits::Rect{2, <: Real})
+    return Makie.apply_transform(axis.transform_func[], Rect2d(limits))
 end
 
+
+
 function axis_setup!(axis::GeoAxis)
     # initialize either with user limits, or pick defaults based on scales
     # so that we don't immediately error
-    targetlimits = Observable{Rect2d}(geodefaultlimits(axis.limits[], axis.source_projection[], axis.target_projection[]))
+    targetlimits = Observable{Rect2d}(geodefaultlimits(axis, axis.limits[]))
     finallimits = Observable{Rect2d}(targetlimits[]; ignore_equal_values=true)
-    transformedlimits = Observable{Rect2d}(finallimits[]; ignore_equal_values=true)
+    inputlimits = Observable{Union{Nothing, Rect2d}}(; ignore_equal_values=true)
     setfield!(axis, :targetlimits, targetlimits)
     setfield!(axis, :finallimits, finallimits)
-    setfield!(axis, :transformedlimits, transformedlimits)
-
-    axis_transform = create_transform(axis.source_projection, axis.target_projection)
+    setfield!(axis, :inputlimits, inputlimits)
 
-    onany(finallimits, axis_transform) do finallims, tfunc
-        transformedlimits[] = Makie.apply_transform(tfunc, finallims)
+    onany(finallimits, axis.transform_func) do finallims, tfunc
+        result = Makie.apply_transform(Makie.inverse_transform(tfunc), finallims)
+        inputlimits[] = result
     end
-    notify(axis_transform)
+    notify(finallimits)
 
     # set up transformed limits
     # TODO get correct values rather than defaulting to all-Earth
 
     topscene = axis.blockscene
-    scenearea = Makie.sceneareanode!(axis.layoutobservables.computedbbox, transformedlimits, DataAspect())
+    scenearea = Makie.sceneareanode!(axis.layoutobservables.computedbbox, finallimits, DataAspect())
     scene = Scene(topscene, px_area=scenearea)
     axis.scene = scene
-    onany(Makie.update_axis_camera, camera(scene), scene.transformation.transform_func, transformedlimits, axis.xreversed, axis.yreversed)
+    onany(Makie.update_axis_camera, camera(scene), scene.transformation.transform_func, finallimits, axis.xreversed, axis.yreversed)
     notify(axis.layoutobservables.suggestedbbox)
     Makie.register_events!(axis, scene)
     on(axis.limits) do mlims

src/utils.jl

@@ -28,7 +28,7 @@ function Makie.apply_transform(t::Proj.Transformation, pt::V) where V <: VecType
     catch e
         # catch this annoying edge case
         # if pt[2] ≈ 90.0f0 || pt[2] ≈ -90.0f0
-        #     println("Caught a 90-lat")
+        #     println("Caught a 90° latitude")
         #     return Point(t(Vec(pt[1], 90.0f0)) ./ PROJ_RESCALE_FACTOR)
         # end
         println("Invalid point for transformation: $(pt)")
@@ -51,10 +51,40 @@ function Makie.apply_transform(f::Proj.Transformation, r::Rect2{T}) where {T}
     end
 end
 
+_apply_inverse(itrans, p) = Makie.apply_transform(itrans, p .* PROJ_RESCALE_FACTOR) .* PROJ_RESCALE_FACTOR
+
 function Makie.inverse_transform(trans::Proj.Transformation)
     itrans = Base.inv(trans)
-    return Makie.PointTrans{2}() do p
-        return Makie.apply_transform(itrans, p) .* PROJ_RESCALE_FACTOR
+    return Makie.PointTrans{2}(Base.Fix1(_apply_inverse, itrans))
+end
+
+function Makie.apply_transform(t::Makie.PointTrans{2, Base.Fix1{typeof(GeoMakie._apply_inverse), Proj.Transformation}}, r::Rect2{T}) where T
+    f = t.f.x
+    xmin, ymin = minimum(r) .* PROJ_RESCALE_FACTOR
+    xmax, ymax = maximum(r) .* PROJ_RESCALE_FACTOR
+    try
+
+        (umin, umax), (vmin, vmax) = Proj.bounds(f, (xmin,xmax), (ymin,ymax))
+
+    if isinf(umin)
+        umin = -180.0
+    end
+    if isinf(umax)
+        umax = 180.0
+    end
+
+    if isinf(vmin)
+        vmin = -90.0
+    end
+    if isinf(vmax)
+        vmax = 90.0
+    end
+
+    return Rect(Vec2(T(umin), T(vmin)),
+                Vec2(T(umax-umin), T(vmax-vmin)))
+    catch e
+        @show r
+        rethrow(e)
     end
 end