Why omit +units?
Closed this issue · 7 comments
Missing +units is also not a bug, you're converting coordinates to meters by using geokeysToProj4.convertCoordinates().
Why did you decide on this format? Probably not the best idea to have a proj4 string that only works with this code
Hello! GeoTIFF's CRS is specified by geokeys, not by proj4. Geokeys are just references to the EPSG database entries, in other words, CRS "parts". EPSG basically allows following cases:
- CRS can use different units for different axes. Although, for now there's no CRS like that, I think it's better be safe than sorry.
- CRS can specify non-linear units. Most notable case is angular units. Angular units may not be specified in degrees, thus conversion is needed. Moreover, CRS can use units such as linear speed and rotation speed for reasons unknown to me. I don't think that any actual image would use such CRS, yet I want to support it.
While I'm not sure if point 2 is handled by +to_meters, point 1 is definitely not.
I'm gonna close this issue, but if you or somebody else have different opinion, feel free to reopen it. I'll also pin it for future references.
i'm not sure if you are aware, but +vunits and +vto_meter also exists in the proj4 spec.
i've also noticed that your +a and +b parameters arent always in meters? you seem to generate them based on the horizontal units. I notice that when i was using these following keys. +a ends up as 6378137/3.28084, instead of just being 6378137
GTCitationGeoKey: NAD83 / WV South (ftUS), GTModelTypeGeoKey: 1, GTRasterTypeGeoKey: 2, GeogAngularUnitsGeoKey: 9102, GeogAzimuthUnitsGeoKey: 9102, GeogCitationGeoKey: NAD83, GeogEllipsoidGeoKey: 32767, GeogGeodeticDatumGeoKey: 32767, GeogInvFlatteningGeoKey: 298.257222101, GeogLinearUnitsGeoKey: 9003, GeogPrimeMeridianGeoKey: 32767, GeogPrimeMeridianLongGeoKey: 0, GeogSemiMajorAxisGeoKey: 6378137, GeographicTypeGeoKey: 32767, PCSCitationGeoKey: NAD83 / WV South (ftUS), ProjCoordTransGeoKey: 8, ProjFalseOriginEastingGeoKey: 1968500, ProjFalseOriginLatGeoKey: 37, ProjFalseOriginLongGeoKey: -81, ProjFalseOriginNorthingGeoKey: 0, ProjLinearUnitsGeoKey: 9003, ProjStdParallel1GeoKey: 38.88333333333334, ProjStdParallel2GeoKey: 37.483333333333334, ProjectedCSTypeGeoKey: 32767, ProjectionGeoKey: 32767, VerticalUnitsGeoKey: 9003
https://proj.org/en/9.3/usage/projections.html#units
... Note that this does not affect the units of linear parameters such as +x_0, +y_0, +a or +b which should always be specified in meters.
i'm not sure if you are aware, but +vunits and +vto_meter also exists in the proj4 spec.
Those are for Z axis. For some reason, I skipped over whole vertical CS thing, I need to refresh my memory.
i've also noticed that your +a and +b parameters arent always in meters?
They should be in meters, looks like a bug.
I'll check these issues at the end of this week, thanks for reporting them!
Hello again! Spec says that GeogSemiMajorAxisGeoKey (and other axes-defining keys) should use units defined in GeogLinearUnitsGeoKey:
In your case units are US feet which can be converted to meters by multiplying by 0.30480060960121924. So +a and +b values are correct, and the problem lies in your file itself.
I'll write back later (today or at the end of the week) about your other question.
P.S. I definitely need to add new aliases, I wasn't aware of newer standard.
Hello! I think I solved Vertical CS thing. GeoTIFF can have following 3D CRS types:
- Geographic 3D - basically [lon, lat] + height.
- Cartesian 3D - [x, y] + height.
- Cartesian 2D + Vertical 1D - [x, y] + height but different.
- Geocentric - [x, y, z] coordinates where Earth's center is at [0, 0, 0] point.
In types 1-3, height (3rd coordinate in proj4) is a band (pixel) value. It's not specified which band to use, but I believe every sane DEM provider will just give you one band.
In types 1-2, we just need to multiply pixel value by the units value which is defined by Cartesian CS. Also, I suppose, presence of keys that define vertical CS is invalid here, thus such keys must be ignored.
In type 3, we need a reference point for the height that we'll just hardcode. Then, I think, we just need to add height of that point to the ellipsoid's axes. And we'll use pixel value multiplied by the units as a Z coordinate.
As for geocentric CRS, I believe, we can just take pixel value as a Z coordinate. I didn't find any examples or other info, so let's think by analogy :)
All of this needs to be implemented, I'm not sure when I'll be done.
In your case units are US feet which can be converted to meters by multiplying by 0.30480060960121924. So +a and +b values are correct, and the problem lies in your file itself.
Ah, I wouldnt doubt it. I'm not actually using an image, i'm using a Lidar file that just so happens to use the same geotiff key format. Definitely possible that somebody wrote the wrong values
Hello! I finally added Z coordinates support, new 2024.03.09 release is on the way. Please note that convertCoordinates() now accepts 4 arguments.