Why omit +units?

[matafokka]

Why omit +units?

[elunty]

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

[matafokka]

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:

1. 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.
2. 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.

[elunty]

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.

[matafokka]

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!

[matafokka]

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.

[matafokka]

Hello! I think I solved Vertical CS thing. GeoTIFF can have following 3D CRS types:

1. Geographic 3D - basically [lon, lat] + height.
2. Cartesian 3D - [x, y] + height.
3. Cartesian 2D + Vertical 1D - [x, y] + height but different.
4. 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.