Feature #396

grdpaste: can't paste geographc grid one with longitudes in [-180,180] range and other in [0,360] range

Added by Eduardo almost 4 years ago. Updated over 3 years ago.

Status:ClosedStart date:2013-10-06
Priority:NormalDue date:
Assignee:-% Done:

0%

Category:-
Target version:Candidate for next minor release
Platform:Linux

Description

I'm making a mosaic of SRTM tiles (from JPL). Since it is a coastal region I must make the missing tiles.

If I make it with

gmt grdmath -R303/304/-36/-35 -I3s -fg 0 = S36W057.hgt

grdpaste says:

grdpaste: Grids do not share a common edge!

but, if I make it with

gmt grdmath -R-57/-56/-36/-35 -I3s -fg 0 = S36W057.hgt

grpaste works fine.

What confuses me first is the info from grdinfo

gmt grdinfo S37W057.hgt
S37W057.hgt: Title: S37W057
S37W057.hgt: Command:
S37W057.hgt: Remark:
S37W057.hgt: Gridline node registration used [Geographic grid]
S37W057.hgt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
S37W057.hgt: x_min: 303 x_max: 304 x_inc: 0.000833333333333 name: x nx: 1201
S37W057.hgt: y_min: -37 y_max: -36 y_inc: 0.000833333333333 name: y ny: 1201
S37W057.hgt: z_min: NaN z_max: NaN name: z
S37W057.hgt: scale_factor: 1 add_offset: 0

gmt grdinfo S36W057.hgt
S36W057.hgt: Title: Produced by grdmath
S36W057.hgt: Command: grdmath -R303/304/-36/-35 -I3s -fg 0 = S36W057.hgt
S36W057.hgt: Remark:
S36W057.hgt: Gridline node registration used [Geographic grid]
S36W057.hgt: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
S36W057.hgt: x_min: 303 x_max: 304 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 1201
S36W057.hgt: y_min: -36 y_max: -35 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 1201
S36W057.hgt: z_min: 0 z_max: 0 name: z
S36W057.hgt: scale_factor: 1 add_offset: 0
S36W057.hgt: format: netCDF-4 chunk_size: 134,134 shuffle: off deflation_level: 0

After a while I've figured out the problem, but I think that grdpaste should be a little more smart.

Associated revisions

Revision 12291
Added by Paul almost 4 years ago

Better longitude checking as in issue #396

History

#1 Updated by Paul almost 4 years ago

  • Status changed from New to Feedback

Eduardo, you did not show what grdpaste command you used, and you did not tell us what you figured out, and in what way grdpaste could be smarter. Was all that left as an exercise for the reader?

#2 Updated by Joaquim almost 4 years ago

OK, let's make a deal and say that grdpaste is half smart

grdmath -R303/304/-36/-35 -I3s -fg 0 = S36W057.grd
grdmath -R-56/-55/-36/-35 -I3s -fg 0 = S36W056.grd
grdpaste S36W057.grd S36W056.grd -Glixo.grd

grdmath -R303/304/-36/-35 -I3s 0 = S36W057.grd
grdmath -R-56/-55/-36/-35 -I3s 0 = S36W056.grd
grdpaste S36W057.grd S36W056.grd -Glixo.grd
grdpaste: Grids do not share a common edge!

but now I'm puzzled because I thought that -fg was of no use in this case.

#3 Updated by Paul almost 4 years ago

The key thing is to use -fg in grdpaste if it was not used in creating the files. The case above works if you add -fg to grdpaste. Without it there is not way for it to know that your files are geographic.

#4 Updated by Joaquim almost 4 years ago

Hmm, right. But I see future others falling in this same trap.
Now I have to find why grdinfo crashes when run on the pasted grid.

#5 Updated by Eduardo almost 4 years ago

I did this:

grdinfo S39W058.hgt
S39W058.hgt: Title: S39W058
S39W058.hgt: Command:
S39W058.hgt: Remark:
S39W058.hgt: Gridline node registration used [Geographic grid]
S39W058.hgt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
S39W058.hgt: x_min: 302 x_max: 303 x_inc: 0.000833333333333 name: x nx: 1201
S39W058.hgt: y_min: -39 y_max: -38 y_inc: 0.000833333333333 name: y ny: 1201
S39W058.hgt: z_min: NaN z_max: NaN name: z
S39W058.hgt: scale_factor: 1 add_offset: 0

(grid at http://sylt.fcaglp.unlp.edu.ar/~esuarez/clase/srtm/)

so I made S40W058.hgt

gmt grdmath -R302/303/-40/-39 -I3s -fg 0 = S40W058.hgt
^^^^^^^^^^^^^
same region

and then:

gmt grdpaste S39W058.hgt S40W058.hgt -GS39_40W058.nc -fg
grdpaste: Grids do not share a common edge!

After a time asking me what I'm doing wrong I've tried

gmt grdmath -R-58/-57/-40/-39 -I3s -fg 0 = S40W058.hgt

and

gmt grdpaste S39W058.hgt S40W058.hgt -GS39_40W058.nc -fg

without error.

#6 Updated by Paul almost 4 years ago

  • Status changed from Feedback to Closed

OK, expected behavior. I will add a word of caution to the grdpaste man page but otherwise this is "user error"

#7 Updated by Eduardo almost 4 years ago

I've read the modification in grdpaste man page:

"Note: For geographical grids, you may have to use -f to handle periodic longitudes unless the input grids are properly recognized as such via their meta-data."

My test grids are recognized as geographic grids (see below), but grdpaste keeps failing.

gmt grdpaste S39W058.hgt S40W058.hgt -GS39_40W058.nc
grdpaste: Grids do not share a common edge!

Where is my error?

================== grdinfo information ======================

grdinfo S39W058.hgt
S39W058.hgt: Title: S39W058
S39W058.hgt: Command:
S39W058.hgt: Remark:
S39W058.hgt: Gridline node registration used [Geographic grid]
S39W058.hgt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
S39W058.hgt: x_min: 302 x_max: 303 x_inc: 0.000833333333333 name: x nx: 1201
S39W058.hgt: y_min: -39 y_max: -38 y_inc: 0.000833333333333 name: y ny: 1201
S39W058.hgt: z_min: NaN z_max: NaN name: z
S39W058.hgt: scale_factor: 1 add_offset: 0

grdinfo S40W058.hgt
S40W058.hgt: Title: Produced by grdmath
S40W058.hgt: Command: grdmath -R302/303/-40/-39 -I3s -fg 0 = S40W058.hgt
S40W058.hgt: Remark:
S40W058.hgt: Gridline node registration used [Geographic grid]
S40W058.hgt: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
S40W058.hgt: x_min: 302 x_max: 303 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 1201
S40W058.hgt: y_min: -40 y_max: -39 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 1201
S40W058.hgt: z_min: 0 z_max: 0 name: z
S40W058.hgt: scale_factor: 1 add_offset: 0
S40W058.hgt: format: netCDF-4 chunk_size: 134,134 shuffle: off deflation_level: 0

#8 Updated by Paul almost 4 years ago

  • Status changed from Closed to In Progress

If you look at the grdinfo report from a geographic grid it will say

S40W058.hgt: Gridline node registration used [Geographic grid]
S40W058.hgt: x_min: 302 x_max: 303 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 1201
S40W058.hgt: y_min: -40 y_max: -39 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 1201

Yours is reported as a geographic grid in one place and Cartesian (x, y) in another. Not sure how you made these:

S39W058.hgt: Gridline node registration used [Geographic grid]
S39W058.hgt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
S39W058.hgt: x_min: 302 x_max: 303 x_inc: 0.000833333333333 name: x nx: 1201
S39W058.hgt: y_min: -39 y_max: -38 y_inc: 0.000833333333333 name: y ny: 1201

Perhaps this is because you are using the ESRI ascii format for the file? It would not have set the units properly.

#9 Updated by Eduardo almost 4 years ago

ha

The problem is in the SRTM file or in grdinfo

grd2xyz S39W058.hgt -fg | minmax
<Standard Input>: N = 1442401 ←58/-57> ←39/-38> ←4/180>

#10 Updated by Paul almost 4 years ago

How did you create the srtm tiles? Did you use the esri format? If so there is nothing in the file that says it is geographic so you would need -fg here.

#11 Updated by Eduardo almost 4 years ago

#12 Updated by Eduardo almost 4 years ago

The grids are in ESRI format, but grdinfo keeps saying "Geographic Grid". Why if it isn't really recognized as geograpic and requires -fg?

Adding -fg gives the same output

gmt grdpaste S39W058.hgt S40W058.hgt -GS39_40W058.nc -fg
grdpaste: Grids do not share a common edge!

#13 Updated by Paul almost 4 years ago

I see now; this is a new feature Joaquim implemented to handle those *.hgt files directly. Because it is a special check, some items in the header are not set correctly so the rest of GMT gets confused. We will see what is needed.

#14 Updated by Paul almost 4 years ago

  • Status changed from In Progress to Resolved

OK, try now, in r12287.

#15 Updated by Eduardo almost 4 years ago

The problem is still there in r12288

#16 Updated by Paul almost 4 years ago

Well, I downloaded the grid from the site you showed me and now grdinfo works for me:

grdinfo S37W057.hgt -M
S37W057.hgt: Title: S37W057
S37W057.hgt: Command:
S37W057.hgt: Remark:
S37W057.hgt: Gridline node registration used [Geographic grid]
S37W057.hgt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
S37W057.hgt: x_min: 303 x_max: 304 x_inc: 0.000833333333333 name: longitude [degrees_east] nx: 1201
S37W057.hgt: y_min: -37 y_max: -36 y_inc: 0.000833333333333 name: latitude [degrees_north] ny: 1201
S37W057.hgt: z_min: -13 at x = 303.079166667 y = -36.3633333333 z_max: 29 at x = 303.225 y = -36.9783333333
S37W057.hgt: scale_factor: 1 add_offset: 0
S37W057.hgt: 9 nodes (0.0%) set to NaN

and with
grdmath -R-58/-57/-37/-36 -I3s -fg 0 = S38W057.hgt

I can do

grdpaste S38W057.hgt S37W057.hgt -Gt.nc

#17 Updated by Eduardo almost 4 years ago

yes! you're doing the same that I did yesterday night!

The only question I have is, how do you know that the correct region is -R-58/-57/-37/-36 if grdinfo reports that the other grid is -R303/304/-37/-36 and no further information.

I think one would tend to use -R302/303/-37/-36 instinctively, if is misguided for that information.

Perhaps you can use the Remark: section of grdinfo to warn users that the .hgt files are in [-180,180] range

Thank you for your patience!

#18 Updated by Paul almost 4 years ago

I think the point is that it does not matter if you use 303/304 or -57/-56 since once -fg is used in grdmath and the hgt file is now understood to be geo then grdpaste knows what to do.

#19 Updated by Eduardo almost 4 years ago

Your example works because the common edge is the latitude.

If the common edge is the longitude grdpaste fails.

Try this:

grdmath -R303/304/-36/-35 -I3s -fg 0 = grid.hgt
grdpaste grid.hgt S37W057.hgt -Gt.nc

#20 Updated by Paul almost 4 years ago

Aahhhh....Will investigate!

#21 Updated by Paul almost 4 years ago

OK, see if this works for you; in r12291.

#22 Updated by Eduardo almost 4 years ago

It works! Thank you!

#23 Updated by Florian over 3 years ago

  • Target version changed from 5.1.0 to Candidate for next minor release

#24 Updated by Remko over 3 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF