Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Why omit +units? #6

Closed
elunty opened this issue Feb 5, 2024 · 7 comments
Closed

Why omit +units? #6

elunty opened this issue Feb 5, 2024 · 7 comments

Comments

@elunty
Copy link

elunty commented Feb 5, 2024

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

@matafokka
Copy link
Owner

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.

@matafokka matafokka pinned this issue Feb 6, 2024
@elunty
Copy link
Author

elunty commented Feb 6, 2024

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
Copy link
Owner

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
Copy link
Owner

matafokka commented Feb 8, 2024

Hello again! Spec says that GeogSemiMajorAxisGeoKey (and other axes-defining keys) should use units defined in GeogLinearUnitsGeoKey:

image

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
Copy link
Owner

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.

@elunty
Copy link
Author

elunty commented Feb 13, 2024

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
Copy link
Owner

Hello! I finally added Z coordinates support, new 2024.03.09 release is on the way. Please note that convertCoordinates() now accepts 4 arguments.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants