tracking projections that work and do not work

[lazarusA]

I did the following subset of projections (almost all available from Proj4) to test how many are done correctly and show them in BeautifulMakie(output). Is better than expected, still some corner cases [that I don't use], but those are really hard to get right. Either way, I leave the code here as reference. This is just for lines, I can imagine that for surfaces (heatmaps) less will be successful.

using GLMakie, GeoMakie
GLMakie.activate!()
let
    projections = ["+proj=adams_hemi", "+proj=adams_ws1", "+proj=adams_ws2",
        "+proj=aea +lat_1=29.5 +lat_2=42.5", "+proj=aeqd", "+proj=airy", "+proj=aitoff",
        "+proj=apian", "+proj=august", "+proj=bacon", "+proj=bertin1953", "+proj=bipc +ns",
        "+proj=boggs", "+proj=bonne +lat_1=10", "+proj=cass", "+proj=cea",
        "+proj=chamb +lat_1=10 +lon_1=30 +lon_2=40", "+proj=collg", "+proj=comill",
        "+proj=crast", "+proj=denoy", "+proj=eck1", "+proj=eck2", "+proj=eck3",
        "+proj=eck4", "+proj=eck5", "+proj=eck6", "+proj=eqc", "+proj=eqdc +lat_1=55 +lat_2=60",
        "+proj=eqearth", "+proj=euler +lat_1=67 +lat_2=75", "+proj=fahey", "+proj=fouc", "+proj=fouc_s",
        "+proj=gall", "+proj=geos +h=35785831.0 +lon_0=-60 +sweep=y", "+proj=gins8", "+proj=gn_sinu +m=2 +n=3",
        "+proj=goode", "+proj=guyou", "+proj=hammer", "+proj=hatano",
        "+proj=igh", "+proj=igh_o +lon_0=-160", "+proj=imw_p +lat_1=30 +lat_2=-40", "+proj=isea",
        "+proj=kav5", "+proj=kav7", "+proj=laea", "+proj=lagrng", "+proj=larr", "+proj=lask",
        "+proj=lcca +lat_0=35", "+proj=leac", "+proj=loxim",
        "+proj=lsat +ellps=GRS80 +lat_1=-60 +lat_2=60 +lsat=2 +path=2", "+proj=mbt_s", "+proj=mbt_fps",
        "+proj=mbtfpp", "+proj=mbtfpq", "+proj=mbtfps", "+proj=merc", "+proj=mill", "+proj=misrsom +path=1",
        "+proj=moll", "+proj=murd1 +lat_1=30 +lat_2=50",
        "+proj=murd3 +lat_1=30 +lat_2=50", "+proj=natearth", "+proj=natearth2",
        "+proj=nell", "+proj=nell_h", "+proj=nicol",
        "+proj=ob_tran +o_proj=mill +o_lon_p=40 +o_lat_p=50 +lon_0=60", "+proj=ocea", "+proj=oea +m=1 +n=2",
        "+proj=omerc +lat_1=45 +lat_2=55", "+proj=ortel", "+proj=ortho", "+proj=patterson", "+proj=poly",
        "+proj=putp1", "+proj=putp2", "+proj=putp3", "+proj=putp3p", "+proj=putp4p", "+proj=putp5",
        "+proj=putp5p", "+proj=putp6", "+proj=putp6p", "+proj=qua_aut", "+proj=robin", "+proj=rouss",
        "+proj=rpoly", "+proj=sinu", "+proj=times", "+proj=tissot +lat_1=60 +lat_2=65", "+proj=tmerc",
        "+proj=tobmerc", "+proj=tpeqd +lat_1=60 +lat_2=65", "+proj=urm5 +n=0.9 +alpha=2 +q=4",
        "+proj=urmfps +n=0.5", "+proj=vandg", "+proj=vandg2", "+proj=vandg3", "+proj=vandg4",
        "+proj=vitk1 +lat_1=45 +lat_2=55", "+proj=wag1", "+proj=wag2", "+proj=wag3", "+proj=wag4",
        "+proj=wag5", "+proj=wag6", "+proj=wag7", "+proj=webmerc +datum=WGS84", "+proj=weren",
        "+proj=wink1", "+proj=wink2", "+proj=wintri"]

    fig = Figure(resolution=(1200, 9000))
    k = 2
    for i in 1:39, j in 1:3
        try
            ga = GeoAxis(fig[i, j]; dest=projections[k],
                title=projections[k], coastlines=true)
            hidedecorations!(ga)
        catch ex
        end
        k += 1
    end
    fig
end;

[asinghvi17]

Neat! This will definitely be useful, I'd love to add it as a test if you're OK with that!

Many of the projections which fail here have singularities or infinities at the traditional lon in (-180,180) and lat in (-90,90) limits. E.g., web mercator cannot show the poles, airy also has singularity conditions. You could try to set lonlims=automatic, latlims=automatic in GeoAxis, which I'm pretty sure should help with at least the Mercator-type and Airy projections.

I'm currently changing the projection system which breaks a lot of this locally (basically all of it, actually) but the only way we could do this in the way which Cartopy does would be if we had knowledge of the domain of any arbitrary projection, which I don't believe Proj gives.

In general, I'd have to say that if lines work, so should surface (for the same domain). Note that all the axes are either well-formed or nonexistent (which happens when the bounds of the axis project to infinity).

[lazarusA]

I'd love to add it as a test if you're OK with that!

Sure, that's why I'm sharing :D

Many of the projections which fail here have singularities or infinities.

Indeed, in this regard maybe some proper initial default values for certain projections and domains will be helpful.

I'm currently changing the projection system which breaks a lot of this locally...

Great!, I hope the user syntax is still similar, which as of now feels natural (in the Makie sense :D ).

[asinghvi17]

The syntax shouldn't change, only the internals :D (right now there are issues with the inverse transformation which are pretty fundamental, for example).

[asinghvi17]

Note for this issue: at this point, with the PR Block geoaxis, there's something strange with the size. This is what I get (notice that the sizes are quite small)

The errors I get are:

Errors

```julia Invalid point for transformation: Float32[-100.0, -5.31714f7] +proj=cass Proj.PROJError("cass: Invalid latitude") ┌ Warning: No strict ticks found └ @ PlotUtils ~/.julia/packages/PlotUtils/bZEEj/src/ticks.jl:191 +proj=cea InexactError(:trunc, Int64, NaN) Invalid point for transformation: Float32[-100.0, -165.28813] +proj=eqdc +lat_1=55 +lat_2=60 Proj.PROJError("eqdc: Invalid latitude") Invalid point for transformation: Float32[-100.0, -164.70897] +proj=euler +lat_1=67 +lat_2=75 Proj.PROJError("euler: Invalid latitude") +proj=geos +h=35785831.0 +lon_0=-60 +sweep=y InexactError(:trunc, Int64, NaN32) +proj=guyou ErrorException("Invalid text boundingbox") +proj=igh ArgumentError("range step cannot be zero") +proj=imw_p +lat_1=30 +lat_2=-40 InexactError(:trunc, Int64, NaN32) Invalid point for transformation: Float32[-100.0, -121.2566] +proj=lcca +lat_0=35 Proj.PROJError("lcca: Invalid latitude") Invalid point for transformation: Float32[-100.0, -171.52461] +proj=murd1 +lat_1=30 +lat_2=50 Proj.PROJError("murd1: Invalid latitude") Invalid point for transformation: Float32[-100.0, -171.27933] +proj=murd3 +lat_1=30 +lat_2=50 Proj.PROJError("murd3: Invalid latitude") +proj=ortho ErrorException("Invalid text boundingbox") Invalid point for transformation: Float32[-100.0, -526.6071] +proj=poly Proj.PROJError("poly: Invalid latitude") Invalid point for transformation: Float32[-100.0, -51990.34] +proj=rouss Proj.PROJError("rouss: Invalid latitude") Invalid point for transformation: Float32[-100.0, -174.89333] +proj=tissot +lat_1=60 +lat_2=65 Proj.PROJError("tissot: Invalid latitude") Invalid point for transformation: Float32[-100.0, -167.77003] +proj=vitk1 +lat_1=45 +lat_2=55 Proj.PROJError("vitk1: Invalid latitude") ```

which leads me to believe that we should: (a) just return NaN or Inf when Proj errors with invalid lat/long. (b) create a custom "text boundingbox finder" which ignores all infinite textboxes, i.e., textboxes with origin at infinity.

[joa-quim]

You may find this PROJ issue informative on why some projections have not worked at all.

[asinghvi17]

Thanks! I encountered that a while ago - that's why the transformation code in utils.jl is so verbose, it's actually a try-catch on each transformation in case Proj reports an error :(.

Maybe we could speed that up in Proj.jl itself, perhaps a fast path when applying a transform (or fast_apply(transform, ...) which returns NaN when the point is invalid for whatever reason?

With the latest changes to geoaxis that I've made, pretty much every projection just works (although we still have quite a bit of work to do to figure out valid limits etc).

[joa-quim]

Sorry, can't comment on that. I know nothing about GeoMakie code. Just thought the info could be useful as I found the same problem when calling PROJ trough GDAL.