What to do if highest resolution grid has dumb increments?

[PaulWessel]

In testing the mars_relief recipe on the highest resolution MOLA grid (Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif) one learns that it is pixel registered and 106694 x 53347 with a 200 m pixel which works out to be 12.1468873601 arc seconds. Not exactly a good number to divide into one degree. Basically, there are 296.372222222 pixels per degree. Hence the trusted tiler which tries to make 10x10 degree JP2 tiles run into massive

grdconvert [WARNING]: (e - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.

warnings and it is just junk of course since one tile cannot match the edge of the next. I had tentatively named this grid mars_relief_12s_p.nc, knowing it is not actually 12s. Of course, anyone studying Mars might want to make the highest quality map they can of a specific region and want that grid, but we are unable to make it tiled. So, the options I see are:

1. Let the highest resolution of a data set be the one with a grid increment that divides into 1 degree and yields a whole integer.
2. Let the highest resolution of a data set be the next standard increment (1,3,5,10,15, ...)
3. Let the tiler script simply skip grids that cannot be tiled and we only serve it as a single grid (as we do for 6m and coarser).

In case 1, we find the first integer below 296.372222222 that divides nicely into 3600 is 288. Thus, one would select 12.5s as the increment and filter the 2.1468873601s grid marginally to yield a 12.5 grid. We may choose to name it to the nearest integer second for compliance with our patterns (so 12s or 13s). The original highest resolution grid would not be distributed.
In case 2, we know the answer is 15s so we simply produce solutions from 15s and up. The original highest resolution grid would not be distributed.
In case 3 we upload the untiled mars_relief_12s_p.nc grid (like we do for low resolutions) and start tiling at 15s. This means anyone attempting to cut a chunk off @mars_relief_12s will have to wait for the entire 3.1 Gb grid to be downloaded (once), but at least the highest resolution grid is distributed.

So the casual user might be fine with cases 1 or 2, while Marsophobes will complain that we messed up the high resolution data by filtering.

I dont like to dumb down the original so I think we should pursue option 3. It is a simple test to add to the tiler to check if we have an integer number of nodes per degree, and if not we skip tiling that guy. In support of this, our unadulterated highest res netCDF grid is 3.1 Gb while the original TIF from NASA is 11 Gb, all due to our lossless compression and use of 16-bit integers with granularity of 50 cm. Comments, please.

[PaulWessel]

PS. I should add that one could probably make tiles that fit but they would not be 10x10 but somewhat variable to ensure we have exact common edges. However, this would mean some work in determining which tiles to get I think as well as the naming of the tiles (030W40N etc would not actually be true)_. But perhaps a bit more thought should go into how to tile something like this as clearly tiling is a superiors solution than to distribute 3.1 Gb grids.

[PaulWessel]

Third post: So 106694 is 2 * 7 * 7621 where the latter is a prime number. Hence one option is to make tiles that are 7621 x 7621 which is roughly 25:42:51.4285713221 degrees on the side. This would then give us 14 x 7 = 98 tiles so now we would have smaller and manageable tiles assuming JP2 is OK with an odd-numbered width. Also, the internal machinery that determines which tiles overlap a region would need to use 25:42:51.4285713221 as the step and have some scheme in translating to tile names (I am guessing we would not call them S38.5714285715E025.7142857143..mars_relief_12s_g.jp but perhaps S38E025..mars_relief_12s_g.jp

[PaulWessel]

Know I am talking to myself here but looks like this could work:

The longitudes of the lower left corner could be reduced to integers this way (we skip the last one at 180):

gmt math -T0/106694/7621 -o1 T 360 106694 DIV MUL FLOOR 180 SUB =
-180
-155
-129
-103
-78
-52
-26
0
25
51
77
102
128
154

while the latitude bands would be

gmt math -T0/53347/7621 -o1 T 360 106694 DIV MUL RINT 90 SUB =
-90
-64
-39
-13
13
39
64

So we would build S39E025..mars_relief_12s_g.jp and while neither 39S, 25E or 12s are correct we have meta data that knows the truth and hence can do the right thing. This way of computing those w/e/s/n values also works in the benign cases of proper increments of course.

[joa-quim]

Know I am talking to myself here

Sorry, it's the day 😄

[PaulWessel]

My 98 tiles seems to work OK - just needed some changes to the tiler script. I think with each tile in the 15-25 Mb range the use of this huge data set will be simple.

[WalterHFSmith]

I think the solution Paul seems to be converging to looks like a good one.
The question of filenames not having coordinates that parse to the exact edge coordinates could be handled by having an optional look-up table so that tiles could be labeled tile_N for N 0 or 1 to whatever, and the look-up table would give the exact coordinates.
I did not like the original suggestion that the originally supplied grid could be "filtered" just a little to generate a grid with nicer step sizes. The reason, not appreciated by most of us who don't know much more than the rudiments of the Nyquist-Shannon sampling theorem, is that even such a "minor" filtering has surprisingly bad impacts on the spectrum of the data, and to surprisingly long wavelengths (long compared to sample step size). I showed this in a paper with Karen Marks where we compared the 2-arc-minute pixel-registered grid of Sandwell & Smith with the ETOPO2 created by others who interpolated the original data onto a 2-arc-minute grid-registered grid. The losses are quite severe at quite long wavelengths.
So I recommend doing something that does not require filtering or interpolating the data, if possible.

Going forward, one can urge people that make data products to consider the issues raised here, and choose a sampling step (if possible) so that (360 degrees) / (sample step X) and (180 degrees) / (sample step Y) are highly composite numbers.

[PaulWessel]

Thanks, I talked myself into the right solution once I realised my script could handle this with a few minor changes.

[Esteban82]

BTW, @PaulWessel now that you are deeling with this.

How is it possible to reuse the machinery (to select the resolution of the remote data set) for datasets on a local pc (see forum for details).

[PaulWessel]

Now that @Esteban82 and I are actually playing with Pluto, I learn that there is no backwards compatible way to handle this since the gmt_data_server.txt entry for the highest resolution Pluto tiles says (and needs to be) 52.0732883317s. Being clever, I had cleverly only assigned 8 bytes for the increment string since I assumed it would be 04s and 10m etc. The overflow deletes the unit and we are left with a crude increment assumed to be in degrees.

The fix in master is trivial: Use 16 or 32 bytes for inc in the structure and all works well (we tested this). Then 6.5 will be required even if you dont care about planets since gmt_data_server.txt gets updated and the parsing above crashes gmt.

To have a solution that works with GMT 6.4, the options without code changes are slim. I think the only thing that is simple is to call gmt_data_server_v3.txt and keep just the parsable ones in gmt_data_server.txt . Then 6.5 will read the _v3 file an 6.x x < 5 will read the original name. The two files only differs by the high-res xx.fffffffs entries. Only 6.5 can access them. Then, going forward we only update v3.txt