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?
- Most of the rest of the questions below are trying to expand
on the nuances of this primary question...
- (I provide more background about my application near the bottom,
but hopefully it is mostly not needed?)
Based on recommendations in documentation, for new
files/imports/etc I plan to switch to WKT1 strings instead
of proj strings, but that doesn't help with previously-imported
compatibility.
QUESTION: What are the potential difference(s) between
tranformations created via different strategies? (Other
question's answers may influence this answer, and/or this answer may
render some of the other questions irrelevant.)
* Plain proj_create() supplied with the from-old-GDAL proj4 string
as-is?
- Currently usually radians, barring pipelines with unit
conversions.
- Like the "Version 4 to 5 API Migration" example.
- One archived email suggests this also might change to degrees
eventually in some future version?
https://lists.osgeo.org/pipermail/proj/2019-April/008427.html
- QUESTION: Are there corner cases that this will not work
as expected? Perhaps non-WGS84 datums?
- See also my speculation about very old pj_fwd()/pj_inv()
approach in application background section below.
* Sequence: proj_create() with black box old string,
proj_crs_get_geodetic_crs() from that,
and proj_create_crs_to_crs_from_pj() from the two?
- Currently usually degrees, I think.
- Like the second of the two quick start for C examples?
- Does this work without "+type=crs"? Documentation suggests
not (not a CRS).
* proj_create_crs_to_crs() with both:
- Currently usually degrees, I think.
- Like the "Version 4 to 6 API Migration" example.
- This seems most similar to what I was doing with PROJ.4,
but abuses(?) the arguments with non-CRSs, and does angular
units differently.
- Presumably: Use a hard-coded argument (something like
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
except maybe with "+type=crs", and/or "something" to do
units conversion steps.)
- And the from-old-GDAL proj4 string?
* Start with from=projCreate("EPSG:4326") or similar,
maybe update it with proj_normalize_for_visualization()
and/or proj_crs_alter_cs_angular_unit(),
and finally combine it with something like
proj_create_ctx_to_ctx_with_pj(ctx,from,proj_create(oldProj4)).
- (Or more likely maybe proj_normalize_for_visualization()
should be applied last to the combination, since it normalizes
both input and output.)
- QUESTION: Given that
https://proj.org/en/stable/operations/operations_computation.html
includes mentions of PROJ looking up special case conversions
between pairs of CRS's in a database, would that be confused
by tweaking the CRS (axes and/or units) before combining
it with the old proj4 string?
- See additional questions about how proj4 strings don't
seem to "actually" be CRS's, units, etc. below.
QUESTION (angular units): What is the "best" way to get radians
for geographic lat/lon coordinates, similar to earlier versions
of PROJ?
* QUESTION: A variation of this question also applies to new
"based on WKT string" option I'm working on. I still want
radians for interacting with existing infrastructure in this
case.
* "Best": Ideally avoid converting back and forth repeatedly,
and also hopefully "simple" to code.
* I originally combined a hard-coded proj4 PJ
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
with the GDAL-supplied proj4 PJ in pj_transform(),
which then usually used radians.
- Now with the new API proj_create_crs_to_crs(), it
usually expects degrees.
* Difficulty: There's a special nuance for cases when
old GDAL OGRSpatialReference::IsGeographic() returned true.
- In these cases, my GDAL plugin would adjust image's
affine (aka world or 6 parameter) linear transform for
radians, since the original affine usually output degrees
in this case.
- QUESTION: What's the bsst way to detect "IsGeographic()"
to account for how that was compensated by the affine part,
ideally without GDAL?
- MAYBE: However, depending on overall strategy above, maybe
I won't need to dig into this?
- MAYBE copy and simplify GDAL IsGeographic() function?
But it seems a tad involved.
- QUESTION: What is a "PJ_TYPE_BOUND_CRS", which seems
relevant to how the new "IsGeographic()" function is
implemented? The documentation lists this enum option,
but doesn't describe it at all. Or is it irrelevant if
the proj4 string came from old GDAL?
* MAYBE just convert angular units myself outside of PROJ?
- But converting back and forth definitely seems a bit
inefficient and ugly.
- And I need to know exactly when to do such conversion...
- QUESTION: Under this option, should I just assume proj
is using degrees? It looks like there are ways it could
(sometimes) be doing something else, although unusual.
- QUESTION: Or should use functions like
proj_angular_input(), proj_angular_output(),
proj_degree_input(), and/or proj_degree_output() to determine
if I should deal with it externally? Any gotchas with
this? What if a CRS is using really unusual angular units
(neither degrees nor radians)?
* MAYBE use proj_create() with the persisted old proj4 string,
instead of proj_create_crs_to_crs()?
- See also question about different ways to create PJ object,
above.
* MAYBE somehow add a units conversion step (pipeline) to the
hard-coded "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
I previously used with pj_transform()?
- See more "projinfo"-based analysis below...
* MAYBE add on an extra units conversion step "later".
- Vaguely like proj_normalize_for_visualization(), except
for units instead of axis order.
- QUESTION: Is there any robust way compose transforms
on objects (instead of blindly playing modification games with
black-box strings, or fully understanding and reproducing
big chunks of the string parser)?
- MAYBE: It looks like proj_crs_alter_cs_angular_unit() sort of
does this.
- Only works with a CRS, not a constructed transform.
- See also above mention, especially the question about
looking up best conversions in the database.
QUESTION: What does proj_create_crs_to_crs() do with non-CRS proj4
args (and similar "expects crs" functions if given a non-CRS)?
* The documentation and examples says it accepts proj4 strings
(without a "+type=crs" qualifier), but doesn't say what it does
with them in this case.
* Testing with the "Version 4 to 6 API Migration example" suggests
non-CRS "works" in the sense it does not error out (returns
non-NULL, no warnings), but testing with "projinfo" suggests
the PJ transform is in a state that projinfo can't tell me
anything about. See also below.
QUESTION: Clarify how additional parameters affect "+proj=latlon"
(especially "+ellps" and "+datum").
* https://proj4.org/en/stable/operations/conversions/latlon.html
* The documentation is good about "on it's own", with:
"No parameters will affect the output of the operation if used
on it's own."
* But I'm not sure what "used in a declaritive manner" means in
"However, the parameters below can be used in a declarative
manner when used with cs2cs or in a transformation pipeline."
- It seems like latlon's "+datum" might alter results if
combined with a projection with a different "+datum" via the
proj_create_crs_to_crs() function?
------------------
Things I have been looking into:
------------------
I'm trying to play with "projinfo" to attempt to learn how things
work. Attempting to summarize (although it is hard):
* In the single proj4 argument case:
- If doing a pipeline or similar non-CRS cases, the WKT
basically just seems to wrap the proj4 args (with
"CONVERSION" and "METHOD", rather than "really" being
native WKT).
- With "+proj=latlon", if I try to use pipeline to convert
input from rad to deg, and/or output from deg to rad, it
errors out with "Mismatched units between step ...".
With or without "+type=crs", and even though in a CRS
context (or two-arg form below) "+proj=latlong" would
appear to use degrees (WKT).
- If I add "no-op" unit conversions to the pipeline
(like "+step +proj=unitconvert +xy_in=rad +xy_out=rad"),
it doesn't seem to optimize them out (the non-CRS thing above).
* In the "-s/-t" two-CRS form, if either or both proj4 strings
lack "+type=CRS" (which they will with old GDAL, at least):
- "projinfo" errors out complaining that it isn't a CRS.
- But proj_create_crs_to_crs() seems to do something
(return non-NULL) even though they aren't CRS's. But
it isn't documented what it does with this non-CRS form,
nor does "projinfo" help.
- QUESTION: Would appending "+type=CRS" to an old GDAL proj4
string be generically valid in all cases? (I would prefer
to treat proj4 as a black box, avoid modifying the
string "blindly", and any manipulation is "as an object" instead
of "as a string".)
- Hard-coded
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +type=CRS"
works. But I haven't been able to enhance that to use a
pipeline to do units conversion.
- QUESTION: I haven't figured out a syntax to make
"+type=CRS" work in a pipeline. Is it possible?
- I've tried starting with the EPSG:4326 WKT string as
one of the arguments, except I modify its CS ANGLEUNITS.
- Depending on the main projection arg, two-CRS "projinfo"
shows a two step process in WKT, but a
simplified "main" single-arg proj4 string. Maybe this
is reasonable? Or does multi-setp WKIT imply it is
slower?
- QUESTION: What if the black-box-from-old-GDAL string happens
to be IsGeographic()? One way to handle that (if
detected - see above) is to convert both input AND output
to radians within PROJ, somehow. (Although there may be
other ways with tweaking old affine...)
- I notice that output from the following
mentions rough "ballback" transformations.
- projinfo -s EPSG:4326 \
-t "+proj=merc +ellps=WGS84 +lat_ts=33 +type=crs"
- Similar for other variations of the 4 to 6 migration
example, including trying to mix its' clrk66 ellipsoid
with WGS84. (I'm not even sure which combinations make
any sense.)
- Adding "+datum=WGS84" to the proj4 string seems to help.
But a "+datum=clrk66" isn't recognized.
- See also documentation bugs in this example below
(mixing up datum vs ellps)?
- I haven't checked: Does old GDAL always (or usually)
include a datum in proj4 strings?
- There is similar confusion trying to mix
---
I tried searching old github issues and mailing list archives, although
it isn't easy, I've been far from thorough, and even if there are
interesting nuggets in there, they are buried under a lot of stuff.
But a few I stumbled over:
https://github.com/OSGeo/PROJ/issues/3546
- Seemed to result in the quick start example
involving proj_crs_get_geodetic_crs().
- https://github.com/OSGeo/PROJ/pull/3515
- https://github.com/OSGeo/PROJ/issues/3508
https://github.com/OSGeo/PROJ/issues/4283
- An angular units issue in cct command line.
https://github.com/OSGeo/PROJ/issues/3597
- Unexpected unit conversion in pipeline, workaround to
explicitly convert to desired units, and apparently proj
has an "optimizer" to eliminate redundant steps.
https://github.com/OSGeo/PROJ/issues/2027
- Querying input/output angular units. (But only degrees and radians,
not arbitrary.)
https://github.com/OSGeo/PROJ/issues/1662
- Seems somewhat related to my questions, but the answers
don't seem great.
https://github.com/OSGeo/PROJ/issues/1178
- Early discussion about axis swapping issues.
- Also mentions units issues, but doesn't really go anywhere
with it.
https://lists.osgeo.org/pipermail/proj/2019-April/008414.html
- The first message in a thread asking about angular units.
----------
DOCUMENTATION BUGS:
----------
Tangent: Things that looked like bugs to me (not just "incomplete")
when reading the documentation:
- BUG https://proj.org/en/stable/operations/conversions/unitconvert.html
Typo: "z_in" is described as "Vertical output units", not input.
- BUG https://proj.org/en/stable/operations/projections/eqc.html
EPSG:1028 and EPSG:1029 are mentioned on this page, but projinfo
doesn't recognize them.
- BUG https://proj.org/en/stable/operations/transformations/deformation.html
The sentence "Both grids GDAL both reads and writes both file formats"
probably needs some editing.
- BUGS: https://proj.org/en/stable/operations/transformations/helmert.html
and various others: Various pages will (sometimes) report
"Erroneous nesting of equation structures" instead of showing a
mathematical expression. (Unless I turn off javascript. But I
don't know TeX or its partial javascript port(s) well enough to
debug these without a lot more effort.)
- BUGS: Examples for version 4 migrations need to include
stdio.h (scanf/printf) and math.h (HUGE_VAL). Also, main()'s
return type should be specified as int rather than implied.
- BUG (or at least confusing): In the "Version 4 to 6 Migration"
section, the original proj4 example program uses a different
proj string than the PROJ6 version of the same program.
- Specifically, datum is changed to ellpse in
"+proj=merc +datum=clrk66 +lat_ts=33" vs
"+proj=merc +ellps=clrk66 +lat_ts=33".
- PROJ 9.8.1 doesn't like the "+datum=clrk66". PROJ 4.9.3's "proj"
program doesn't seem to like it either. Has it ever been valid?
- QUESTION: If this is actually "correct", then how would a
programmer know when/how to change the string like this?
----------
APPLICATION BACKGROUND:
----------
20 years ago I wrote some code for handling background image and
projection code on the Earth for a proprietary application.
The application draws the Earth as a WGS-84 shaped ellipsoid, with
heights of vertices adjusted based on terrain, and images texture
mapped onto it. Then it draws other things in the same coordinate
system, using Earth to provide context.
- The application has "native" support for this that doesn't
involve GDAL or PROJ at all. Most of this support thinks in
either ECEF or (lat,lon,alt) in *radians* and meters.
- GDAL and PROJ are isolated in plugins and custom abstractions
to support a lot more "standard" image file formats and
the projections specified in those files, assuming the
plugins are available.
- Known bug: Currently it treats "MSL altitude" from various sources
as exactly the same thing as "height above the ellipsoid", without
handling even a low-precision geoid (vertical datum) as appropriate.
Roughly up to a 100 meter vertical error in some cases.
This has been open in our issue tracker for nearly a decade, and
I've known about it for far longer. (Despite the decade delay,
this is probably higher priority than a couple of meters from
things like continental drift.)
The application includes an ability to "import" imagery into a
application-native format, with generated pyramid / RRDs to
handle any eye point / zoom level.
- I'm glossing over lots of image stack composition,
multiresolution pyramid / RRD data structures, and eyepoint logic.
- Currently the projection in such an imported image is based
on proj4 strings as generated by GDAL ca. 2.3.1
OGRSpatialReference::exportToProj4(), plus the
affine (aka world or 6 parameter) transform that most ground
images use.
- Actually, I adjust the affine values slightly very early
for "normalized" resolution-independent image coordinates
(range 0 to 1 for opposite sides of the image) and
sometimes radians vs degrees difficulty (next bullet).
- Slight difficulty: My GDAL plugin had some additional
logic to test OGRSpatialReference::IsGeographic(), and
if true adjust the returned affine transform values for
degrees vs radians, since in this case the original affine
values are generally referenced against degrees instead of
radians.
- To convert between geographic lat/lon and (normalized) image
coordinates, the plugin used pj_transform() between the
"black box" proj4 string originally from GDAL, and a hard-coded
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
along with my own code for the affine transform (or its
inverse). (ECEF conversion happens elsewhere.)
- During development 20 years ago, I initially used pj_fwd()
and pj_inv(), but within a few days changed to pj_transform()
under the unconfirmed speculative theory that pj_transform()
might have a chance of handling datum changes?
- QUESTION: Is there any truth to this theory? I haven't seen
hard evidence either way...
---------
Matthew Ogilvie
---------
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
