[gdal-dev] +proj=wintri +over option works weel in gdalwarp but not in ogr2ogr

[Joaquim Manuel Freire Luís]

Eve, thanks for looking to this but I think I’ve not been clear enough in my goals.
I don’t want to clip into the [0 360] interval. Quite on the contrary, I want to make maps that extend that limit but keep ofc the earth periodicity. My initial attempts worked quite well with the Van Der Grinten projection (https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/82?u=joaquim).

For the Winkel Triple gdalwarp actually works quite well too

try this with the grid in the issue https://github.com/OSGeo/gdal/issues/7644

gdalwarp -t_srs "+proj=wintri +over" earth15.grd lixo.tiff

It’s the vector step that fails (but the cl540 has 384 segments, not just one).

Anyway, this is not an important matter. It would just be nice if I could to with “+proj=wintry +over” what I can with “+proj=vandg +over”

Note, I saw the warning at https://proj.org/usage/projections.html#longitude-wrapping

“Note however that for most projections where the 180 meridian does not project to a straight line, +over will have no effect or will not lead to expected results.”

The Van Der Grinten doesn’t fulfill that condition but things still work for it, so I don’t really know how to interpret that warning.

Joaquim

From: Even Rouault <even.rouault at spatialys.com>
Sent: Sunday, April 30, 2023 2:43 PM
To: Joaquim Manuel Freire Luís <jluis at ualg.pt>; gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] +proj=wintri +over option works weel in gdalwarp but not in ogr2ogr


Joaquim,

I've looked at that. Experimentally, it doesn't seem that the current maths for wintri in +over mode give reasonable results for longitudes > 360deg in absolute value (and I'm not sure if there would be a natural extension possible. Far beyond my knowledge of the underlying maths).

gdalwarp does a little better by default because it couples the forward and inverse transformer to check the domain of validity and experimentally must find this limitation and will thus give up for source pixels whose longitude > 360 deg. For ogr2ogr / vector this is more complicated as reprojection of a single geometry  (here cl540.gpkg is a single polygon with longitudes in [-540,540] range). is a "all points succeed to reproject => OK" / "one single point fails to reproject => KO" type of situation, at least in the logic followed by OGRGeometry::transform() logic.

So while the following patch to PROJ is probably technically correct, it wouldn't help in practice

diff --git a/src/projections/aitoff.cpp b/src/projections/aitoff.cpp
index 9d060999f..cce319157 100644
--- a/src/projections/aitoff.cpp
+++ b/src/projections/aitoff.cpp
@@ -59,6 +59,11 @@ static PJ_XY aitoff_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
     struct pj_opaque *Q = static_cast<struct pj_opaque *>(P->opaque);
     double c, d;

+    if( fabs(lp.lam) > 2 * M_PI ) {
+        proj_errno_set(
+                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
+        return xy;
+    }
     c = 0.5 * lp.lam;
     d = acos(cos(lp.phi) * cos(c));
     if (d != 0.0) { /* basic Aitoff */

The only solution I can offer is that you must explictly clip your source geometry to the validity range of wintri with -clipsrc -360 -90 360 90.

Even


Le 30/04/2023 à 03:11, Joaquim Manuel Freire Luís via gdal-dev a écrit :
Even, this is a kind of continuation of the subject that I brought up in https://github.com/OSGeo/gdal/issues/7644

The file size limitation of attachments doesn’t help to make this case easily reproducible, but the story is, I can make that rectangular map using gdalwarp and ogr2ogr that use the +over option and +proj=vandg
(see https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/82?u=joaquim)

But when I try the same with +proj=wintry the gdalwarp op works well but the ogr2ogr doesn’t.
(see https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/85?u=joaquim)

You can more less reproduce this case using this file
http://fct-gmt.ualg.pt/tmp/cl540.gpkg

ogr2ogr  -t_srs "+proj=wintri +over" cl540_wintri.gpkg cl540.gpkg

Why would that be in this case?

Joaquim



_______________________________________________

gdal-dev mailing list

gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230430/fabd8148/attachment-0001.htm>