On Thu, May 21, 2026 at 09:53:39AM -0600, Matthew Ogilvie wrote:
> PRIMARY QUESTION: Given a previously imported/converted image file
> with a proj4 string from GDAL ca. version 2.3.1 or earlier
> embedded in it, what is the "best" way to create an appropriately
> compatible transformation object with new PROJ versions like 9.8.1?
In the absense of any responses, I'm tentatively using (roughly)
the following code snippet.
Maybe someone will alert me if anything seems non-ideal in it,
or if good, maybe it could be an example to help anyone else with
similar questions.
========= CUT ===========
PJ *wgs84=proj_create(ctx,"EPSG:4326");
wgs84=proj_normalize_for_visualization(ctx,wgs84);
wgs84=proj_crs_alter_cs_angular_units(ctx,wgs84,"radian",1.0,NULL,NULL);
////////
PG *newCrs;
if(oldProj4)
{
/* Adapted from pj_add_type_crs_if_needed() internally used
* by proj_create_crs_to_crs() only (not other forms of creation).
*/
std::string adjStr=proj4String;
if( ( starts_with(adjStr, "proj=") ||
starts_with(adjStr, "+proj=") ||
starts_with(adjStr, "+init=") ||
starts_with(adjStr, "+title=") ) &&
adjStr.find("type=crs") == std::string::npos )
{
adjStr.append(" +type=crs");
}
newCrs = proj_create(ctx,adjStr.c_str());
PJ_TYPE pjType=proj_get_type(newCrs);
if( pjType==PJ_TYPE_GEOGRAPHIC_2D_CRS ||
pjType==PJ_TYPE_GEOGRAPHIC_3D_CRS )
{
newCrs = proj_crs_alter_cs_angular_unit(ctx,newCrs,
"radian",1.0,NULL,NULL );
}
}
else // new WKT
{
newCrs = proj_create(ctx,wktString);
}
newCrs = proj_normalize_for_visualization(ctx,newCrs);
////////
PJ *transform;
transform = proj_create_crs_to_crs_from_pj(ctx,wgs84,newCrs,NULL,NULL);
// Then just use that final "transform"...
========= CUT ===========
Notes:
- Obviously the snippet above is leaving out a lot of error checking,
memory leak cleanup, and other tangential details. (Real code has
those...)
- Most piecemeal construction / tweaking of a transform apparently
requires actual CRS's, which needs the " +type=crs" appended to the
proj4 string, or else later CRS-dependent steps fail.
1. Except for plain proj_create_crs_to_crs(), which does that for
you. Perhaps this could be clarified in the documentation?
2. MINOR CONCERN: As mentioned in my original email, I don't think
this is very robust (I don't think it works with pipelines,
etc), but it is probably good enough for my purposes here.
- CONCERN: Are there more PJ_TYPE's I need to worry about?
- I notice that GDAL's OGRSpatialReference::IsGeographic() may
also do a bunch more stuff for "bound CRS"s that I don't
really understand.
- FYI: If I unconditionally do newCrs's
proj_crs_alter_cs_angular_unit() adjustment and print out
the proj4String and WKT representations of the final
resulting transform, it looks like non-geopraphic
newCrs's will create transforms that involve extra steps
that aren't optimized properly. (As near as I can tell.)
- The documentation mentions that finding the best CRS to CRS
transform sometimes involves looking up special cases in
the proj database file.
- https://proj.org/en/stable/operations/operations_computation.html
- CONCERN: Could customizing axes (order, units) mess
this up in some cases when it does such lookups in the final
proj_crs_to_crs_with_pj() call?
- ALTERNATIVE: In most (all?) cases it might work to
just use "transform = proj_create(proj4String);" directly (no
explicit crs_to_crs at all).
- CONCERN: Could this ever result in a different transform?
Lack of datum conversions, etc? If different, how should
I decide which would be better?
- Future CONCERN: One of the archived emails I linked to from
my first email mentioned the idea of changing this to work
in degrees as well?
- FYI: My original email mentions these concerns and others
in more detail.
- FYI: My original email also has a section about what looks
like documentation bugs / typos someone may want to fix...
---------
Matthew Ogilvie
---------
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
