Bug #541

grdpaste: changes in dX and/or dY after paste

Added by Eduardo about 3 years ago. Updated almost 3 years ago.

Status:ClosedStart date:2014-04-16
Priority:NormalDue date:
Assignee:-% Done:

90%

Category:-
Target version:Candidate for next bugfix release
Affected version:5.x-svn Platform:Linux

Description

I was playing with CGIAR's SRTM tiles. I know, its format is awfull.

But...

grdinfo srtm_22_19.asc=ei -fg
srtm_22_19.asc: Title:
srtm_22_19.asc: Command:
srtm_22_19.asc: Remark:
srtm_22_19.asc: Pixel node registration used [Geographic grid]
srtm_22_19.asc: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
srtm_22_19.asc: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
srtm_22_19.asc: y_min: -35.0004161824 y_max: -29.9995828491 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 6001
srtm_22_19.asc: z_min: NaN z_max: NaN name: z
srtm_22_19.asc: scale_factor: 1 add_offset: 0

grdinfo srtm_22_20.asc=ei -fg
srtm_22_20.asc: Title:
srtm_22_20.asc: Command:
srtm_22_20.asc: Remark:
srtm_22_20.asc: Pixel node registration used [Geographic grid]
srtm_22_20.asc: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
srtm_22_20.asc: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
srtm_22_20.asc: y_min: -40.0004163035 y_max: -34.9995829701 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 6001
srtm_22_20.asc: z_min: NaN z_max: NaN name: z
srtm_22_20.asc: scale_factor: 1 add_offset: 0

grdpaste srtm_22_19.asc=ei srtm_22_20.asc=ei srtm_22_1920.grd -fg

grdinfo srtm_22_1920.grd
srtm_22_1920.grd: Title: srtm_22_1920.grd
srtm_22_1920.grd: Command: grdpaste srtm_22_19.asc=ei srtm_22_20.asc=ei -Gsrtm_22_1920.grd -fg
srtm_22_1920.grd: Remark:
srtm_22_1920.grd: Pixel node registration used [Geographic grid]
srtm_22_1920.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
srtm_22_1920.grd: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
srtm_22_1920.grd: y_min: -40.0004163035 y_max: -29.9995828491 y_inc: 0.000833333343422 name: latitude [degrees_north] ny: 12001
srtm_22_1920.grd: z_min: -22 z_max: 6928 name: z
srtm_22_1920.grd: scale_factor: 1 add_offset: 0
srtm_22_1920.grd: format: netCDF-4 chunk_size: 131,130 shuffle: on deflation_level: 3

Note how dY changed from 0.000833333333333 to 0.000833333343422 to accommodate discrepancy.

But the things are even worse

grdinfo srtm_22_21.asc=ei -fg
srtm_22_21.asc: Title:
srtm_22_21.asc: Command:
srtm_22_21.asc: Remark:
srtm_22_21.asc: Pixel node registration used [Geographic grid]
srtm_22_21.asc: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
srtm_22_21.asc: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
srtm_22_21.asc: y_min: -45.0004164245 y_max: -39.9995830912 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 6001
srtm_22_21.asc: z_min: NaN z_max: NaN name: z
srtm_22_21.asc: scale_factor: 1 add_offset: 0

grdpaste srtm_22_1920.grd srtm_22_21.asc=ei -Gsrtm_22_1921.grd -fg

grdpaste says nothing about the different dY, even a warning.

And now dY changed again!

grdinfo srtm_22_1921.grdsrtm_22_1921.grd: Title: srtm_22_1920.grd
srtm_22_1921.grd: Command: grdpaste srtm_22_1920.grd srtm_22_21.asc=ei -Gsrtm_22_1921.grd -fg
srtm_22_1921.grd: Remark:
srtm_22_1921.grd: Pixel node registration used [Geographic grid]
srtm_22_1921.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
srtm_22_1921.grd: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
srtm_22_1921.grd: y_min: -45.0004164245 y_max: -29.9995828491 y_inc: 0.000833333346785 name: latitude [degrees_north] ny: 18001
srtm_22_1921.grd: z_min: -79 z_max: 6928 name: z
srtm_22_1921.grd: scale_factor: 1 add_offset: 0
srtm_22_1921.grd: format: netCDF-4 chunk_size: 131,129 shuffle: on deflation_level: 3

If I repeat the paste with the next tile dY changes again to 0.000833333348467

History

#1 Updated by Joaquim about 3 years ago

Eduardo, I have spent countless hours debugging cases that turned out to be shit coming out from CGIAR tiles that never learned the difference between grid vs pixel registration nor how to compute the tile limits correctly after the grid to pixel (needless) conversion.

I did not check it now but I strongly suspect that the problem is still the same. The grid (pixel reg) limits are not correct (within the possible of the double precision variables) and that small error propagates into DY.

#2 Updated by Eduardo about 3 years ago

Then I suspect that is better not to use CGIAR grids. I'll tell that to my colleagues here.

#3 Updated by Joaquim about 3 years ago

Aside from the gaps filling work that was nicely done, I never understood why one would use a re-interpolated (with the corresponding loss inherent to the process) version of the originals.

#4 Updated by Joaquim about 3 years ago

I wonder if grdpaste should not automatically detect those 1x1 deg pix reg tiles at 3 arc-sec resolution and simply ignore the limits in the header and use the exactly known (their double representation) limits.

#5 Updated by Eduardo about 3 years ago

I don't know, maybe to ban those grids. I did a grdreformat of both ESRI and GeoTiff grids and the results are not the same. How can I tell this is the right grid?
@
grdreformat srtm_22_21.asc=ei grid_a.grd -fg
grdreformat srtm_22_21.tif=gd grid_b.grd -fg

grdinfo grid_a.grd
grid_a.grd: Title: Produced by grdreformat
grid_a.grd: Command: grdreformat srtm_22_21.asc=ei grid_a.grd -fg
grid_a.grd: Remark:
grid_a.grd: Pixel node registration used [Geographic grid]
grid_a.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
grid_a.grd: x_min: 284.999583479 x_max: 290.000416812 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
grid_a.grd: y_min: -45.0004164245 y_max: -39.9995830912 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 6001
grid_a.grd: z_min: -79 z_max: 3332 name: z
grid_a.grd: scale_factor: 1 add_offset: 0
grid_a.grd: format: netCDF-4 chunk_size: 131,131 shuffle: on deflation_level: 3

grdinfo grid_b.grd
grid_b.grd: Title: Produced by grdreformat
grid_b.grd: Command: grdreformat srtm_22_21.tif=gd grid_b.grd -fg
grid_b.grd: Remark:
grid_b.grd: Gridline node registration used [Geographic grid]
grid_b.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
grid_b.grd: x_min: 285.000000145 x_max: 290.000000145 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 6001
grid_b.grd: y_min: -44.9999997579 y_max: -39.9999997579 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 6001
grid_b.grd: z_min: -79 z_max: 3332 name: z
grid_b.grd: scale_factor: 1 add_offset: 0
grid_b.grd: format: netCDF-4 chunk_size: 131,131 shuffle: on deflation_level: 3
@

#6 Updated by Joaquim about 3 years ago

There might be a GMT hand on this. Grids imported via GDAL are always converted to grid line registration due to the absence on information from GDAL on the original registration type. See issue #254

#7 Updated by Eduardo about 3 years ago

Aha, so what is the Best Current Practice in dealing with this grids? What is your advice?

#8 Updated by Joaquim about 3 years ago

At this moment a rock solid solution would be to grdedit them to be grid registered and impose the known limits (easy, they fall on integer numbers). But we need a definitive solution for this.

#9 Updated by Joaquim almost 3 years ago

  • Status changed from New to Resolved

On r13222 and r13226 I implemented the capability of reading and writing the GeoTiff key(s) that inform whether a grid is pixel or grid registered. This means we now can write pixel and grid registered and same for reading. There will be still an ambiguous case for files that have not set the AREA_OR_POINT key. In this case, the old default (grids are always grid reg) is maintained.
As far as I can tell, this the most we can do to address this issue and (#543). So I'm declaring this solved ... till next round

#10 Updated by Joaquim almost 3 years ago

  • Status changed from Resolved to Closed
  • % Done changed from 0 to 90

Also available in: Atom PDF