WGS84 ellipsoid not properly converted
Closed this issue · 12 comments
My geotiff file has these geokeys:
GTModelTypeGeoKey: 1
GTRasterTypeGeoKey: 1
GeogEllipsoidGeoKey: 7030
GeogGeodeticDatumGeoKey: 6326
PCSCitationGeoKey: "Transverse Mercator; WGS84; WGS84"
ProjCoordTransGeoKey: 1
ProjFalseEastingGeoKey: 500000
ProjFalseNorthingGeoKey: -5300000
ProjLinearUnitsGeoKey: 9001
ProjNatOriginLatGeoKey: 0
ProjNatOriginLongGeoKey: 19
ProjScaleAtNatOriginGeoKey: 0.9993
ProjectedCSTypeGeoKey: 32767
ProjectionGeoKey: 32767
When passing it through toProj4, the proj4 string is returned like this:
+proj=tmerc +lat_0=0 +lon_0=19 +k_0=0.9993 +x_0=500000 +y_0=-5300000 +a=6378137 +b=6378137 +no_defs
First of all, this doesn't look right: when +a and +b are equal, doesn't that mean that the ellipsoid is a sphere?
EPSG:7030 is the WGS84 ellipsoid, which is definitely not a sphere.
Second, passing this string to the proj4 library actually throws an exception:
Incorrect elliptical usage. Try using the +approx option in the proj string, or PROJECTION["Fast_Transverse_Mercator"] in the WKT.
I see that the definition being output is from GeogEllipsoidGeoKey.js, but I don't know why that definition is like that.
In EllipsoidsNamesToProj.js there is a definition for ellipsoid 7030, and it is mapped to "WGS84" - so I would have expected the output to be +ellps=WGS84.
In any case, the mapping seems to be wrong:
https://epsg.io/7030-ellipsoid.gml
clearly states an inverse flattening.
For reference, gdalsrsinfo reports the SRS of the geotiff file like this:
PROJ.4 : +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +datum=WGS84 +units=m +no_defs
OGC WKT2:2018 :
PROJCRS["Transverse Mercator; WGS84; WGS84",
BASEGEOGCRS["WGS_1984",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["unknown",0.0174532925199433]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",19,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9993,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",-5300000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Hmm, could it perhaps be that in line 19 of EllipsoidsWorker.js:
let fromCode = names[result.name.toString()];
result.name should be result.id? I don't have a copy of the database, so I'm just guessing based on the iogp-ellipsoid-7030.gml file.
Thank you for reporting this issue and such informative description.
This behavior is indeed wrong. Yes, there might be a mistake in ellipsoids worker, easy to mess something up with such database structure. I also suspect it might be related to the string building process (main.js file).
Anyway, thank you for everything, I appreciate your efforts and interest in this project.
As I said, I'll try to fix it tomorrow. Also, I'll check what's up with that "Incorrect elliptical usage" error, even though it's not related to your main issue.
I just had another look at the GeogEllipsoidGeoKey.js file. Something must have gone wrong when you recreated it, as for all ellipsoids +a and +b are equal!
Goodness, previous release had the same thing! The problem's gotta be in data preparation. How does it even pass my tests?! I'll get more samples just in case stuff like this popped up somewhere else.
Hello, fixed the thing. Your suggestion was right, and there also was a typo at the end of ellipsoids worker which caused your issue. I feel dumb right now.
Please, temporarily switch version in package.json to github:matafokka/geotiff-geokeys-to-proj4#master and see if your samples work. Later today I'll look for more samples and check if there's more same issues.
By the way, can I ask you for your samples? I'll add them to the tests, it surely will help with debugging.
Just pushed a new version to npm, everything seems fine to me. Please, let me know if there's any more bugs.
The 7030 ellipsoid is still not correctly converted. I now get this output:
+proj=tmerc +lat_0=0 +lon_0=19 +k_0=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +a=6378137 +b=6399521.685754821 +no_defs
This cannot be right, since the b term is larger than the a term. b is the minor axis (the distance from the earth center to the poles), so it should be smaller than a (radius at the equator). The correct minor axis for the WGS84 ellipsoid is 6356752.31424518.
The bug is in line 43 of EllipsoidsWorker.js:
b = a / result.f + a;
The correct formula should be
b = a - a / result.f;
because the inverse flattening f is a/(a-b) (https://en.wikipedia.org/wiki/Flattening)
I wonder if it's really a good idea to provide both the +ellps and the +a +b terms. It seems proj relies on +a +b instead of using the predefined +ellps ellipsoid. (As long as the parameters match, that is of course not a serious problem.)
Please, temporarily switch version in
package.jsontogithub:matafokka/geotiff-geokeys-to-proj4#masterand see if your samples work. Later today I'll look for more samples and check if there's more same issues.
Thanks for this cool tip - I didn't know one could directly refer to a github repo from package.json, this is neat!
Thank you, just checked, sign is really off, let me do a quick fix...
wonder if it's really a good idea to provide both the +ellps and the +a +b
I think so. There're ellipsoids with slight modification, for example, Clarke 66 and Clarke 66 Michigan. I'm not sure if +ellps defines anything else rather than +a and +b, so I thought it's good idea to put both and let Proj4 resolve it. Secondly, geokeys might override axes, so you might end up with something like +ellps=... +a=... other_params....
Thanks for this cool tip
No problem :D This also works for a specific version, useful for getting a version that was on GitHub, but wasn't on npm when it was released. Used it to polyfill color input with ColorJS.
Done, and fixed another issue related to the ellipsoids. Now it returns expected value. Hope ellipsoids are finally fine now :D
Now it seems to do the right thing for my files :-)
👍
Awesome! Thank you for using my stuff and for your help, professor!