Horizontal striping in earth_relief_30s?

[Esteban82]

There are horizontal lines in offshore area. I think it is an error in the processing.

This is 15s without those lines.

gmt grdimage @earth_relief_30s -png Earth_30s -I -RFK+r1 -Baf
gmt grdimage @earth_relief_15s -png Earth_15s -I -RFK+r1 -Baf

[maxrjones]

@PaulWessel, do you think this is worthwhile to look at before remaking the earth_relief files?

[PaulWessel]

Yes, we should do that. I suspect the filter width is just a tad too small so that alternate rows ends up with uneven results. Maybe a simple test would be to take a grdcut of the 15s exampe above and just do the grdfilter step that the recipe would call for (i.e., first reproduce the same plot but manuallyl via grdfilter) then increase the filtering a bit. This is what we have

FILTER_WIDTH=$(gmt math -Q ${SRC_RADIUS} 2 MUL PI MUL 360 DIV $INC MUL 0.05 ADD 10 MUL RINT 10 DIV =)

Perhaps the 0.05 is too small - try 0.1. For 30s the above yields

1. Perhaps first try to replace RINT with CEIL and use 0.1. Then we get 1.1 instead of 1.

[PaulWessel]

@Esteban82, I take it there are no problems in the 01m and coarser grids? if so then I think it will go away with another padding than 0.05 and with CEIL rather than RINT.

[Esteban82]

@Esteban82, I take it there are no problems in the 01m and coarser grids?

Exactly. I only find the problem with that 30s resolution.

[PaulWessel]

First check on filtering the original grid to 30s with the same parameters as the script we us gave no stripes. So now running the full job to see if it gives stripes (the files on the server has stripes obviously) since it all looks good on paper. FYI, the 30s grid is filtered with a 5x5 node filter which should be perfect for 15s going to 30s and it looked fine when I do it manually...

[PaulWessel]

Maybe something to do with our rounding to nearest 0.5 meter and/or compression - we will find out.

[PaulWessel]

I can confirm the srv_downsampler_grid.sh does create those stripes; not clear why yet.

[PaulWessel]

It is starting to smell of 64-bit issues. So far this is only reproducible if I filter the global 15s grid and request an output region that is almost the entire globe (the production script does the entire globe). For instance, requesting output for -R-180/180/-90/80 gives the same problem, but -R-180/180/-90/75 does not. But it cannot have to be the input grid size since all the filterings are done on the same grid and only the 30s shows this issue. Of course, holding both the 15s and 30s grids in memory is quite a bit of memory but I have 64 Gb so should not be an issue....

[PaulWessel]

I cannot find a specific problem, but likely to be related to int64 counting. I was able to do a work-around: Filter the N and S hemispheres separately, grdpaste them, then run grdedit to set the remark and then grdconvert to shrink to 50cm ints:

gmt grdfilter SRTM15_V2.4.nc -Fg1 -D1 -I30s -rp -Gnew_30s_p_n.grd--PROJ_ELLIPSOID=Sphere -R-180/180/0/90
gmt grdfilter SRTM15_V2.4.nc -Fg1 -D1 -I30s -rp -Gnew_30s_p_n.grd --PROJ_ELLIPSOID=Sphere -R-180/180/0/90
gmt grdpaste new_30s_p_n.grd new_30s_p_s.grd -Gboth.grd -V
gmt grdedit both.grd -D+r"Reduced by Gaussian Cartesian filtering (1 km fullwidth) from SRTM15_V2.4.nc [Sandwell et al., 2022; https://doi.org/10.1029/2021EA002069]"
gmt grdconvert both.grd -Gearth_relief_30s_p.grd=ns+s0.5+o0 --IO_NC4_DEFLATION_LEVEL=9 --IO_NC4_CHUNK_SIZE=4096

The result does not have any of those stripes when you run

gmt grdimage -R62W/57W/53S/51S earth_relief_30s_p.grd -Cgeo -B -I -png new

@meghanrjones, I think it is best to do this for all the 30s products derived from a 15s product, so probably gebco, gebco_si as well as the Scripps files.

[maxrjones]

Is this something that should be included in the recipes for those datasets?

[PaulWessel]

Hm, well I had hoped to find and fix the bug but it has proven elusive. Perhaps I should add this special trick to the srv_downsampler_grid script. I can do that since we may be doing this for a while if we cannot find the problem,

[PaulWessel]

I have opened a PR for you @meghanrjones. I think I got the region and pieces correctly (the steps above had cut/paste issues).

[maxrjones]

I have opened a PR for you @meghanrjones. I think I got the region and pieces correctly (the steps above had cut/paste issues).

Thanks. I try out remaking the datasets tomorrow.

[PaulWessel]

Just leaving a note here that the recent fix in GMT regarding macOS vDSP use for scaling large output grids does not affect this - it is clearly a grdfilter issue.

[PaulWessel]

Hi @Esteban82, perhaps you can have a look at this and maybe zoom into other areas to see if I am right. Below are two bathymetry plots with no shading, using the area -R57:40W/57W/48:10S/47:50S with the SRTM15 file and then the 30s global filtered file (with the stripes above);

15s original:

30s filtered:

Now, we apply -I+d as default shading:

15s with shading:

30s with shading:

To my eyes, the 15s -> 30s filtering looks reasonable: The pixels are twice as large and hence it looks a bit blockier but the transitions look reasonable and nothing odd stands out. Once we do shading then the 15s looks OK while the 30s looks very striped. This leads me to wonder if the problem lies with the grdgradient call that grdimage does in order to do the shading. What do you think, @GenericMappingTools/core ?

[PaulWessel]

Maybe create some N-S cross-profiles of the 15s and 30s to see if there are any odd steps in 30s.

[Esteban82]

Yes, I was thinking on making profiles.

[Esteban82]

First I made these 2 figures. When I compare the 15s (_p by default) vs 30s_p the profiles (a NS black line in the middle of the map) loks very similar.

However, I got this very strange map when I used 30s_g.

[Esteban82]

gmt grdimage @earth_relief_30s_g -png Earth_30s_g -I -R-57:40/-57/-48:10/-47:50 -Baf -B+t"Gridline"
gmt grdimage @earth_relief_30s_p -png Earth_30s_p -I -R-57:40/-57/-48:10/-47:50 -Baf -B+t"Pixel"

[Esteban82]

This is GEBCO

gmt grdimage @earth_gebco_30s_g -png Gebco_30s_g -I -R-57:40/-57/-48:10/-47:50 -Baf -B+t"Gridline"

[Esteban82]

I reopened by mistake.

[WalterHFSmith]

Hi,

I am glad you were able to close the issue, but I didn’t understand what the resolution was.

Did someone mistakenly compare the GEBCO product to a product from the SRTM15Plus or SYNBATH of Sandwell et al?

Or was the pointer to the 30s grid selecting the wrong grid?

Is there an issue with the way the filtering/downsampling is done at high latitude that makes gradients look strange in cylindrical projections?

Thanks,

W

On Jan 16, 2023, at 9:07 AM, Federico Esteban @.***> wrote:

I reopened by mistake.

— Reply to this email directly, view it on GitHub https://github.com/GenericMappingTools/remote-datasets/issues/32#issuecomment-1384114698, or unsubscribe https://github.com/notifications/unsubscribe-auth/APUT6GMNL2ICLVXFA7F65C3WSVIY3ANCNFSM5QPSYXBQ. You are receiving this because you are on a team that was mentioned.

[Esteban82]

Hi, I am glad you were able to close the issue, but I didn’t understand what the resolution was.

The issue was original closed before we restart the discussion. By mistake, I re-open the issue.

I am not sure if the problem was solved or if I find a new problem. I'm a bit lost and maybe I messed up with the data servers.Therefore I closed the issue as it was previously.-

Could someone reproduce the issue with the 30s_g grids?

[PaulWessel]

gmt grdimage @earth_relief_30s_g -png Earth_30s -I -RFK+r1 -Baf

gives pretty much the same plot at the top of this thread, so that part seems OK.

[Esteban82]

OK, so I don't know what is my problem. I remove the grids and I got a wrong figure. I use oceania server. What server should I use?