Add square conformal projections from libproject

[kbevers]

This PR adds five new projections to PROJ:

adams_hemi: Adams Hemisphere in a Square
adams_ws1: Adams World in a Square I
adams_ws2: Adams World in a Square II
guyou: Guyou
peirce_q: Pierce Quincuncial

The code originates from Gerry Evendens libproject and has been adapted
to work with modern PROJ. To ensure that the modified code works as
intended extensive test data has been created using libproject and
sproj to ensure that no errors occured when porting from libproject to PROJ.
The test data is wrapped in a gie files. All test cases reproduce
results from libproject at the mm level.

Docs are still missing at this point. Will be added before merging.

This work is a first step towards implementing the Spilhaus projection as requested in #1851.

• Closes #xxxx
• Tests added
• Added clear title that can be used to generate release notes
• Fully documented, including updating docs/source/*.rst for new API

[mwtoews]

I see that adams_wsI and adams_wsII are the names that Gerry originally had (look here), however this naming scheme with upper case Roman numerals is different than any other projection name, which are all lowercase with Arabic numerals (e.g. wink1 for "Winkel I", or find them all with git grep "PROJ_HEAD(").

Should these projection names be adams_ws1 and adams_ws2 to follow similar projection names in use?

[kbevers]

Should these projection names be adams_ws1 and adams_ws2 to follow similar projection names in use?

That's a good question. Consistency is probably better than libproject compatibility (don't think it was ever used much).

[nyalldawson]

This is great stuff, nice work @kbevers!

[rouault]

Looking at https://raw.githubusercontent.com/Esri/projection-engine-db-doc/master/csv/pe_list_projcs.csv, one can see the following:

54098,54098,"PE_PCS_WORLD_ADAMS_SQ_2","WGS_1984_Adams_Square_II","PROJCS[""WGS_1984_Adams_Square_II"",GEOGCS[""GCS_WGS_1984"",DATUM[""D_WGS_1984"",SPHEROID[""WGS_1984"",6378137.0,298.257223563]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]],PROJECTION[""Adams_Square_II""],PARAMETER[""False_Easting"",0.0],PARAMETER[""False_Northing"",0.0],PARAMETER[""Scale_Factor"",1.0],PARAMETER[""Azimuth"",0.0],PARAMETER[""Longitude_Of_Center"",0.0],PARAMETER[""Latitude_Of_Center"",0.0],PARAMETER[""XY_Plane_Rotation"",0.0],UNIT[""Meter"",1.0]]","Adams Square II","Esri","10.8.0","no","World",-90.0,90.0,-180.0,180.0
54099,54099,"PE_PCS_WORLD_SPILHAUS_OCEAN_MAP","WGS_1984_Spilhaus_Ocean_Map_in_Square","PROJCS[""WGS_1984_Spilhaus_Ocean_Map_in_Square"",GEOGCS[""GCS_WGS_1984"",DATUM[""D_WGS_1984"",SPHEROID[""WGS_1984"",6378137.0,298.257223563]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]],PROJECTION[""Adams_Square_II""],PARAMETER[""False_Easting"",0.0],PARAMETER[""False_Northing"",0.0],PARAMETER[""Scale_Factor"",1.0],PARAMETER[""Azimuth"",40.17823482],PARAMETER[""Longitude_Of_Center"",66.94970198],PARAMETER[""Latitude_Of_Center"",-49.56371678],PARAMETER[""XY_Plane_Rotation"",45.0],UNIT[""Meter"",1.0]]","Spilhaus World Ocean Map in a Square","Esri","10.8.0","no","World",-90.0,90.0,-180.0,180.0

So Esri has added an explicit XY_Plane_Rotation parameter, that they set to 0 for the WGS_1984_Adams_Square_II, and to 45 for Spilhaus.
libproject hardcodes 45 for . Should we have a rot_xy parameter (like in healpix ?). That parameter actually could potentially be handled in a generic way for all projections.

[busstoptaktik]

I do not have much to contribute, but a few "speeling missteaks" is probably also worth correcting...

[kbevers]

So Esri has added an explicit XY_Plane_Rotation parameter, that they set to 0 for the WGS_1984_Adams_Square_II, and to 45 for Spilhaus.

The Adams square II is actually subject to a 45 rotation as seen here:

PROJ/src/projections/adams.cpp

Lines 174 to 178 in 472712e

	if (Q->mode == ADAMS_HEMI || Q->mode == ADAMS_WSII) { /* rotate by 45deg. */
	double temp = xy.x;
	xy.x = RSQRT2 * (xy.x - xy.y);
	xy.y = RSQRT2 * (temp + xy.y);
	}

So for the Spilhaus perhaps this should just be a 90 degree rotation? Or 0? Not sure of the direction here...

But yeah, a generic rotation parameter would also work. Are you aware of other cases where a parameter like that is used by Esri or others?

[rouault]

So for the Spilhaus perhaps this should just be a 90 degree rotation? Or 0? Not sure of the direction here...

Hum, I'm not sure. A possibility would be that xy_rot=0 would skip the current rotation (if that's what Esri WGS_1984_Adams_Square_II does ?), and xy_rot=45 should apply it (Spilhaus use case ?). Would be good if someone could check against a Esri install (I don't have one) what they get for the above 2 definitions (and possibly modify the base WGS_1984_Adams_Square_II one, to just test different values of the rotation parameter), so we are consistent with what they do.

Are you aware of other cases where a parameter like that is used by Esri or others?

Digging into Esri public repo of definitions, I see that they use it at least for Krovak (where it can be 0 or 90, as long as X_Scale and Y_Scale to apply axis reversal) and Oblique Mercator (where it maps to PROJ +gamma parameter)

[nyalldawson]

Would be good if someone could check against a Esri install

Can you clarify exactly what you're after? I can test tomorrow

[rouault]

Can you clarify exactly what you're after? I can test tomorrow

For example, take a world map (Natural Earth or whatever) and reproject it to ESRI:54098 "WGS_1984_Adams_Square_II" and ESRI:54099 "WGS_1984_Spilhaus_Ocean_Map_in_Square"

And if that's possible test reprojecting to the following custom variations

PROJCS["WGS_1984_Adams_Square_II_with_XY_rotation_20",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Adams_Square_II"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Azimuth",0.0],PARAMETER["Longitude_Of_Center",0.0],PARAMETER["Latitude_Of_Center",0.0],PARAMETER["XY_Plane_Rotation",20.0]

PROJCS["WGS_1984_Adams_Square_II_with_XY_rotation_45",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Adams_Square_II"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Azimuth",0.0],PARAMETER["Longitude_Of_Center",0.0],PARAMETER["Latitude_Of_Center",0.0],PARAMETER["XY_Plane_Rotation",45.0]

[kbevers]

I think this PR is more or less finished now. Thanks for all the input, it's been super helpful!

[nyalldawson]

@rouault

, take a world map (Natural Earth or whatever) and reproject it to ESRI:54098

and ESRI:54099 "WGS_1984_Spilhaus_Ocean_Map_in_Square"

PROJCS["WGS_1984_Adams_Square_II_with_XY_rotation_20",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Adams_Square_II"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Azimuth",0.0],PARAMETER["Longitude_Of_Center",0.0],PARAMETER["Latitude_Of_Center",0.0],PARAMETER["XY_Plane_Rotation",20.0]

It refuses to import either of these custom WKT strings

[rouault]

@nyalldawson OK, so comparing images, at least for ESRI:54098 "WGS_1984_Adams_Square_II", with XY_Plane_Rotation=0, they do apply the 45° final rotation that PROJ also applies. Thus XY_Plane_Rotation must be an extra rotation on top of that. Interesting that it refuses to import custom WKT with other XY_Plane_Rotation values. Smells like only ESRI:54098 and ESRI:54099 cases have been somehow hardcoded.

[nyalldawson]

@rouault there's always a chance I was importing wrong. There's zero documentation for this and no error messages - just no result after the import. ESRI documentation is just plain woeful.