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

Can't determine if CRS is correct #59

Open
sagimann opened this issue Feb 15, 2023 · 3 comments
Open

Can't determine if CRS is correct #59

sagimann opened this issue Feb 15, 2023 · 3 comments

Comments

@sagimann
Copy link

I have cropped a tif file using latest elevation lib, e.g.:

eio clip -o dem.tif --bounds 35.079732 30.901217 35.845468 31.870474

Then inspected the CRS of the file using rasterio, which returned "EPSG:4326", plus some point inside the area (lat,lng), which returned 69:

import rasterio
with rasterio.open('my.tif') as dataset:
    print(f'CRS: {dataset.crs}') # output EPSG:4326
    row, col = dataset.index(lng, lat)
    elevation = dataset.read(1, window=((row, row+1), (col, col+1)))
    print(elevation[0, 0]) # output 69

I'm new to this field, but looking up EPSG:4326 says it is "WGS 84, latitude/longitude coordinate system based on the Earth's center of mass, used by the Global Positioning System among others."

Does this mean the generated tif file contains heights above the WGS84 ellipsoid? I'm confused because we have an independent and highly-accurate measuring system that explicitly reports 69m above MSL and NOT above the ellipsoid. And that system also tells us they are ~20m apart at that point, so they can't both be right..

What am I missing here?

@gdt
Copy link

gdt commented Feb 15, 2023

4326 is a 2-dimensional coordinate system so it does not make sense for elevation of any kind of be expressed in it.

The real question is where dem.tif came from, what the CRS is of that, and what the documentation/metadata says. Clipping should not really change that.

In WGS84 there are two concepts:

  • height above ellipsoid (HAE) which is purely geometric and unrelated to gravity. Unless you are doing surveying or geodesy you probably don't want this
  • WGS84 orthometric height: The distance above the geoid model that is part of WGS84. Today, that means EGM2008, but people often use EMG96 anyway, even though it is formally obsolete. People also use truncated representations of the model, leading to errors in the output (compared to the full model which is by definition correct for WGS84 orthometric height).

@gdt
Copy link

gdt commented Feb 15, 2023

Also, "above MSL" is a fuzzy term and I am surprised if anything "highly accurate" is using it. Around the earth the average sea level varies in height relative to the geoid. Ask them what they really mean. Perhaps it is WGS84 orthometric height -- but it is difficult to measure that precisely. Perhaps they mean ITRF HAE converted via EGM2008, which is almost the same thing. Perhaps they mean some national height datum. Most of those are very close to WGS84 orthometric height, at the meter or two level, and to each other. There is one in Europe (Holland?) which is ~2m different.

I should have said earlier that a DEM that was in HAE, while conceptually sensible, seems very odd to me.

For a much longer treatise on height, see https://opencommons.uconn.edu/cgi/viewcontent.cgi?article=1000&context=nrme_monos&httpsredir=1&referer=

@sagimann
Copy link
Author

sagimann commented Mar 3, 2023

The real question is where dem.tif came from, what the CRS is of that, and what the documentation/metadata says. Clipping should not really change that.

As mentioned, dem.tif is the output of the eio command line provided by this project. If my understanding is correct, it comes from SRTM which provides data in both ellipsoid and EGM96 frames of references. The problem is gdalinfo says the file's CRS is WGS84 but when querying the elevation values, I get EGM96 values instead. So either I don't understand the header value or eio gets the EGM96 variant from NASA but puts a WGS84 header on it...

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