Format - SRTM HGT

[Metatron]

According to the NASA documentation:

3.0 Data Formats The names of individual data tiles refer to the longitude and latitude of the lower-left (southwest) corner of the tile (this follows the DTED convention as opposed to the GTOPO30 standard). For example, the coordinates of the lower-left corner of tile N40W118 are 40 degrees north latitude and 118 degrees west longitude. To be more exact, these coordinates refer to the geometric center of the lower left sample, which in the case of SRTM3 data will be about 90 meters in extent.

Though in the loader it is expected that the coordinate definition lies on the corner of the sample, not on it's center and tries to fix it accordingly. That is wrong - the following would comply with the NASA definition:

GDALDataset* SRTMHGTDataset::Open(GDALOpenInfo* poOpenInfo)

poDS->adfGeoTransform[0] = southWestLon; poDS->adfGeoTransform[1] = 1.0 / (numPixels-1); poDS->adfGeoTransform[2] = 0.0000000000; poDS->adfGeoTransform[3] = southWestLat + 1; poDS->adfGeoTransform[4] = 0.0000000000; poDS->adfGeoTransform[5] = -1.0 / (numPixels-1);

Furthermore:

SRTM1 data are sampled at one arc-second of latitude and longitude and each file contains 3601 lines and 3601 samples. The rows at the north and south edges as well as the columns at he east and west edges of each cell overlap and are identical to the edge rows and columns in the adjacent cell.

Though in the loader is is expected to cover 1.0 degree on the [0,1200]/[0,3600] or [0,1201)/[0,3601) interval. That is wrong, it is slightly more (the overlap). 1.0 degree is then on the [0,1199]/[0,3599] or [0,1200)/[0,3600) interval - the following would correct that also:

poDS->adfGeoTransform[0] = southWestLon; poDS->adfGeoTransform[1] = 1.0 / (numPixels-2); poDS->adfGeoTransform[2] = 0.0000000000; poDS->adfGeoTransform[3] = southWestLat + 1; poDS->adfGeoTransform[4] = 0.0000000000; poDS->adfGeoTransform[5] = -1.0 / (numPixels-2);

I detected the error by trying to connect independent SRTM-tiles in a 3D-software, they weren't adjacent. I implemented the fix, and now the data fits absolutely correct.

Ciao

Niels

[Even Rouault]

The point you're raising is always a difficult question, but I'm not convinced your fix is right and I'm not sure there's really a problem in the existing code.

I've tried your fix though on N43E002.hgt and N43E003.hgt. With your fix, gdalinfo reports :

On N43E002.hgt :
Origin = (2.000000000000000,48.000000000000000)
Pixel Size = (0.000834028356964,-0.000834028356964)
Corner Coordinates:
Upper Left  (   2.0000000,  48.0000000) (  2d 0'0.00"E, 48d 0'0.00"N)
Lower Left  (   2.0000000,  46.9983319) (  2d 0'0.00"E, 46d59'53.99"N)
Upper Right (   3.0016681,  48.0000000) (  3d 0'6.01"E, 48d 0'0.00"N)
Lower Right (   3.0016681,  46.9983319) (  3d 0'6.01"E, 46d59'53.99"N)
Center      (   2.5008340,  47.4991660) (  2d30'3.00"E, 47d29'57.00"N)

On N43E003.hgt :
Origin = (3.000000000000000,48.000000000000000)
Pixel Size = (0.000834028356964,-0.000834028356964)
Corner Coordinates:
Upper Left  (   3.0000000,  48.0000000) (  3d 0'0.00"E, 48d 0'0.00"N)
Lower Left  (   3.0000000,  46.9983319) (  3d 0'0.00"E, 46d59'53.99"N)
Upper Right (   4.0016681,  48.0000000) (  4d 0'6.01"E, 48d 0'0.00"N)
Lower Right (   4.0016681,  46.9983319) (  4d 0'6.01"E, 46d59'53.99"N)
Center      (   3.5008340,  47.4991660) (  3d30'3.00"E, 47d29'57.00"N)

So, you can see that the Lower Left corner is not located on rounded coordinates as you expected. Secondly, instead of having one row/column of overlapping, you have now 2 row/column of overlapping. You can see it with openev easily.

The current geotransform computation makes one row/column of overlapping on each edge, which is expected. One alternative to the current implementation could be too shift the whole grid from half the pixel size to make sure that the south west corner is on rounded coordinates, but I'm not sure it is what must be done.

[Metatron]

Yes I saw the second issue was bogus, Sorry.

I found verification of the first issue by overlaying the SRTM3 v3 GeoTIFFs with raw SRTM v2 HDRs. In SRTM3 pixel-centers are on full Degrees, not pixel-corners.

I don't understand your output, it is saying that the N43E002 tile is on North 48 East 2? Also I think it's not very difficult to understand (if you get the intervals rights, as I didn't because of all the trees :-), NASA says 1 degree, that's it.

Overlaying SRTM3 v2 and v3 reveals the pixel-corner/pixel-center issue.

Ciao

Niels

[Even Rouault]

(Sorry, in my previous post, I messed up the name of my HGTs, they were N47E002.hgt and N47E002.hgt)

I've read again and again the sentences :

For example, the coordinates of the lower-left corner of tile
N40W118 are 40 degrees north latitude and 118 degrees west
longitude. To be more exact, these coordinates refer to the
geometric center of the lower left sample, which in the case
of SRTM3 data will be about 90 meters in extent.

My understanding of this is that : the coordinates of the center of the lower left pixel/sample are N40 W118. So, the lower left extremity of the dataset is a bit (half the pixel size, 0.00041667 degree ~= 46m) southern and western of N40 W118.

Thus, what is reported currently seems correct :

gdalinfo N47E002.hgt:

Origin = (1.999583333333333,48.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (   1.9995833,  48.0004167) (  1d59'58.50"E, 48d 0'1.50"N)
Lower Left  (   1.9995833,  46.9995833) (  1d59'58.50"E, 46d59'58.50"N)
Upper Right (   3.0004167,  48.0004167) (  3d 0'1.50"E, 48d 0'1.50"N)
Lower Right (   3.0004167,  46.9995833) (  3d 0'1.50"E, 46d59'58.50"N)
Center      (   2.5000000,  47.5000000) (  2d30'0.00"E, 47d30'0.00"N)

Furthermore, a bit later the ​spec says :

SRTM3 data are sampled at three arc-seconds and contain
1201 lines and 1201 samples with similar overlapping
rows and columns. This organization also follows the DTED convention.

Well, my understanding is that this means that gdalinfo on a SRTM HGT and the equivalent DTED should report the same coordinates, which is currently the case.

gdalinfo N47E002.dt1:

Origin = (1.999583333333333,48.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (   1.9995833,  48.0004167) (  1d59'58.50"E, 48d 0'1.50"N)
Lower Left  (   1.9995833,  46.9995833) (  1d59'58.50"E, 46d59'58.50"N)
Upper Right (   3.0004167,  48.0004167) (  3d 0'1.50"E, 48d 0'1.50"N)
Lower Right (   3.0004167,  46.9995833) (  3d 0'1.50"E, 46d59'58.50"N)
Center      (   2.5000000,  47.5000000) (  2d30'0.00"E, 47d30'0.00"N)

Maybe you could provide a sample of the different files you mentionned that do not overlay properly so we can sort out this better. Are they supposed to overlay exactly anyway... (maybe some convention has changed between the versions) ?

I'd be interested in the opinion of a "referee"... Frank ?

[Metatron]

My understanding of this is that : the coordinates of the center of the lower left pixel/sample are N40 W118. So, the lower left extremity of the dataset is a bit (half the pixel size, 0.00041667 degree ~= 46m) southern and western of N40 W118.

That is correct - fetching the data at (0,0) has to give me the height of N40W118 and not N39.9..W118.0.. That's the whole issue. You assume "pixels" but there are no pixels, there are height-values, and they lie on a grid defined by the convention.

The SRTM v3 are GeoTIFFs, they are available as 5x5 merged and reworked SRTM v2 tiles. These are my reference files:

SRTM v3:

6000x6000 without overlap, 4.99.93 degree, N15W090 ​http://telascience.sdsc.edu/tela_data/SRTM/srtm.csi.cgiar.org/srtm_v3/SRTM_Data_GeoTiff/Z_19_10.TIF

SRTM v2:

1201x1201 with 1 pixel overlap, 1 degree, N14W090 ​ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM3/South_America/N14W090.hgt.zip

[warmerdam]

The GDAL geotransform is based on an origin at the top left corner of the top left pixel. GDAL *ignores* the idea that some products have pixels without area (ie point-as-pixel) so even for these SRTM datasets the GDAL geotransform origin should be offset by half a pixel from the "location" of the point elevation. That is, the area reported by gdalinfo should be "out half a pixel size" from the apparent region for stuff like the SRTMv2 data.

I'm not clear on differences between v2 and v3 data. It seems from Metatron's last comment that the SRTM v3 data does not have overlapping postings on the boundary so perhaps should be handled differently.

Even - if you would like to take this bug report (reassign to yourself) that would be fine with me! I don't really have time to dig into these SRTM issues deeply.

PS. Isn't the srtmhgt driver new in trunk? I am surprised to see "version: 1.4.2" listed here.

[Even Rouault]

This is a 1.5.0-dev ticket indeed

[Even Rouault]

Metatron,

I've downloaded your 2 files and opened them into OpenEV. Effectively, there's a half pixel size shift between the two. I've also opened them into Global Mapper 8 and it gives the same result as OpenEV. Of course, it's not a formal proof that GDAL is doing the right thing...

Furthermore, you can read the following notice in the ​readme.txt :

PROCESSED SRTM DATA VERSION 3

The data distributed here are in ARC GRID, ARC ASCII and Geotiff format,
in decimal degrees and datum WGS84.  They are derived from the
USGS/NASA SRTM data.  CIAT have processed this data to provide
seamless continuous topography surfaces.  Areas with regions of no
data in the original SRTM data have been filled in using interpolation
methods.  A full technical report on this method is in preparation.

So, SRTMV V3 has clearly been reworked and they maybe decided to follow another convention as the SRTM V2 convention. The whole key of the story is that SRTM V2 has overlapping, where as SRTM V3 has not.

For SRTM V2 (HGT), in GDAL, when you want the elevation at a given location, you find the pixel inside which this location is. For N47E002.hgt, to find the height of N47:00:00 E2:00:00, you must get the pixel whose bottom left coordinates are N46:59:58.5, E1:59:58.5 and upper right coordinates are N47:00:01.5, E2:00:01.5.

For SRTM V3 (GTIFF), in GDAL, when you want the elevation at N47:00:00 E2:00:00, you must get the pixel whose bottom left coordinates are N47:00:00, E1:59:57, and upper right coordinates are N47:00:03, E2:00:00

If you absolutely want the two products to be perfectly overlayed, you have to apply the half pixel size shift But I have the feeling that the SRTM V2 convention is more natural than the SRTM V3 convention.

I've tried your first proposed geotransform :

  poDS->adfGeoTransform[0] = southWestLon; 
  poDS->adfGeoTransform[1] = 1.0 / (numPixels-1);
  poDS->adfGeoTransform[2] = 0.0000000000;
  poDS->adfGeoTransform[3] = southWestLat + 1;
  poDS->adfGeoTransform[4] = 0.0000000000;
  poDS->adfGeoTransform[5] = -1.0 / (numPixels-1);

which gives the perfect overlaying. However as a side effect, it means than we have the column at the east and the line at the south which are completely outside the degree coverage. Nothing in the SRTM V2 spec says that the overlapping between adjacent tiles is not 'symetrical'.

If you still disagree with this analysis, I think you should contact the producers of data/authors of specifications of the both versions to get from first-hand the right way to interpret the data...

[Metatron]

No I don't disagree, I mean the consensus is now clear.

But I see it problematic, that you have to find out the type of data to interprete the values. Both files are height-fields on a well defined grid. Wouldn't it streightforward for GDAL to realize, that elevations are pixel-centered whereas textures are pixel-cornered? (which is without problems, putting the texture-map over the terrain will show no flaw) What for me is even more problematic is that I have to learn and think, that the same value is in two different locations (Greenwich is not in Westend).

But then again: what is the solution for the application? I have to gain knowledge of all these kinds of file-formats and conventions, and apply my own transformations to align them? Isn't that the supposed work of GDAL, taking that away?

Or, after all we assume the SRTM v3 makers made an error, I did not check any other files yet. But in the long run I'm going to mix SDTS, SRTM, GTOPO30 and ArcGrid formats. Let's see how much more trouble is coming.

[Even Rouault]

Currently, GDAL does not try to be clever about the semantics of the data (it doesn't make the difference between an height map, radiometric data, satellite imagery, etc etc). It has a unified model for retrieving pixels (that have an area, not grids), dealing with geotranform matrix and projections. There is not even an API to get a pixel value at a location. You must use the geotransform matrix to compute the pixel you are interest in. Maybe it should evolve.

I think you could launch a thread on gdal-dev list to see what other people think about this issue. Maybe someone will come with an idea / API to deal with that.

I don't know if the SRTM v3 makers made an error. I just tend to think that the convention used in SRTM V2 is easier to use with current GDAL data model because pixels are seen as small areas and not points on a grid.

[warmerdam]

The normal mechanism for indicating that a dataset is "point" samples is to set the GDALMD_AREA_OR_POINT metadata item to GDALMD_AOP_POINT (as opposed to the default area). This does not influence the geotransform which will still effectively treat the raster as pixels with area for georeferencing purposes.

Is this sufficient for the above dilema?

[Even Rouault]

As suggested by Frank and to be consistant with DTED driver that serves the same purpose as SRTMHGT, I've set the GDALMD_AREA_OR_POINT metadata item to GDALMD_AOP_POINT in r18117.

By the way there was a more general debate in ticket #2178 about that whole topic of pixel convention vs area convention