image shading problem for tiled earth relief data?

[seisman]

Plotting a global earth relief map using the following command:

gmt grdimage -JN120/20c -Rg -Baf @earth_relief_05m -I+d -pdf map

It looks good at first glance, but after zooming in, I saw a lot of weird curves:

I don't see these curves for @earth_relief_06m or earth_relief_05m without shading.

[PaulWessel]

Wow, strange. Well ,good that it is not in the data files I guess. Running grdgradient manually and passing the file works. So something with passing @file to grdgradient via GMT_Call_Module.

[PaulWessel]

Just an update. This only applies to tiled data sets. The 06m and larger do not have this problem.

[PaulWessel]

And of course has nothing to do with modern mode. A plain

gmt grdimage -JN120/20c -Rg -Baf @earth_relief_05m -I+d > t.ps

yields the same problem.

[PaulWessel]

Some more experimentation.:

1. I tried a non-curvilinear projection -JQ. No problems, so this involves the projection machinery in some way.
2. I added a gmt_M_memset call to first set the entire intensity grid to zero after projection. The resulting map had the stripes, so they are in the data grid.
3. I then flipped that to instead set the dataset to zero and plotted. The resulting map had intensity (shading) only and it was fine with no stripes.
4. I made a copy of the data grid before passing it to grdgradient, then restored it. Made no difference.

Conclusions:

1. Since the intensity grid is derived from the topo grid via an internal call to grdgradient, and since the intensity grid is fine, then any problem with the data grid must happen during projection.
2. The argument against happening afterwards is that the data and intensity grid take exactly the same route through the projection steps, so it is unclear how the projection could mess up the data grid but not the intensity grid.

Yet, will need to see if there is any difference in the two calls to gmt_grd_project.

[PaulWessel]

Slowly making progress on understanding this but not out of the woods. Can @seisman and @joa-quim try the same command but turning off anti-aliasing? That gives a good result for me but need to verify before I dig into the why (since it works find with a single file or without shading...)

gmt grdimage -JN120/20c -Rg -Baf @earth_relief_05m -I+d -n+a -pdf map

[seisman]

Yes, this command gives me the correct plot.

[PaulWessel]

OK, so I will try to see what is different with the default anti-aliasing when input was tiles versus not, since this has always worked just fine for all cases. Those curved lines has a spacing that is proportional to the grid spacing, are projection dependent, and only appear when auto-shading is on and main file is a tile....

[PaulWessel]

Same error for pixel or gridline tiles, BTW.

[PaulWessel]

Basically, these are the two steps that grid projection takes:

1. Basically do a blockmean: Project each input node, determine which output mode it falls in, and add to sum and number of points. Some output nodes will have many values, others will have 1, some have none.
2. Reverse: Loop over all the projected nodes, recover the nearest lon/lat, and interpolate what the value should be there from the input grid, then either add taht to the sum or use as is.

The -n+a controls this. Default is to do 1 + 2. +a turns off step 1 and only relies on step 2.

The reason we do all this is that step 2 is great if we in effect are interpolating, but would alias if the projected input points are denser than the output points. Step 1 deals with that by averaging (filtering), and then at the end we blend. The result may vary continuously across the map depending on projection, and has worked will for GMT always.
The current problem must introduce some condition we are blind to.

[PaulWessel]

Something related to central meridian and range:

This shows the error (even -JN179/20c shows a little bit):

gmt grdimage -JN0/20c -Baf @earth_relief_05m -I+d -R0/360/80/90 -pdf map05

This one has no issues:

gmt grdimage -JN180/20c -Baf @earth_relief_05m -I+d -R0/360/80/90 -pdf map05

Have a look at this one:

gmt grdimage -JN90/20c -Baf @earth_relief_05m -I+d -R0/360/80/90 -png map05

You can see the lines with land-colored pixels around lon = 60 must be coming from the lat around long = -60. WHy this only happens in bands is the mystery since clearly most pixel end up looking OK.

[seisman]

Looks similar to issue GenericMappingTools/pygmt#515.

[PaulWessel]

Was that with tiled data or full-grid data?

[seisman]

It's 01d full-grid data, but is passed to GMT as an xarray.

[PaulWessel]

I hacked grdimage to write out the number of projected points that fall into each rectangular bin (step 1 above - the blockmean) and then made a plot of it. I did this for tiles and the single 5m pixel grid. The plots are identical so only plosting one here:

Read the cpt as green: 0 hits (outside area), 1, 2, up to 6 max. The pattern reflects the curved paths seen in the final image. However, the final image, which combines the block-meaned points with interpolated points, looks fine for the single grid. So will try to visualize step 2.

[PaulWessel]

Most of the image is blue (1 hit) with the dominant lines from before being red (2 hits). Step 2 code says this:

		z_int = gmt_bcr_get_z (GMT, I, x_proj, y_proj);

if (!GMT->common.n.antialias || nz[ij_out] < 2)	/* Just use the interpolated value */
	O->data[ij_out] = (gmt_grdfloat)z_int;
else if (gmt_M_is_dnan (z_int)) {		/* Take the average of what we accumulated */
	if (nz[ij_out])
		O->data[ij_out] /= nz[ij_out];	/* Plain average */
	else
		O->data[ij_out] = GMT->session.f_NaN;
}
else {					/* Weighted average between blockmean'ed and interpolated values */
	inv_nz = 1.0 / nz[ij_out];
	O->data[ij_out] = (gmt_grdfloat) ((O->data[ij_out] + z_int * inv_nz) / (nz[ij_out] + inv_nz));
}

So, for the blue and red areas this means:

• blue: One hit. Here we simply use the BCR spliine interpolated value z_int
• red: Two hits: Compute wighted average of 3 points: two (the sum with combined weight 2 and z_ing with weight 1/2).

As n goes up the z_in is given less and less weight. E.g., if 5 points inside the sum then the z_int weight is 1/5 vs 5. The idea behind this was to avoid aliasing by favoring the sum if there were many points ( here meaning lots of 5m gridlines going through the same projected node, so splining the lon/lat grid at this point would badly alias the solution), whereas at the other end (no hits) all we do is spline and that is fine. But none of this has changed in years and it has worked fine, except for the tiles. Clearly, the red and blue nodes above give very different results to the point it stands out with the DEM cpt.

[PaulWessel]

Another observation is that towards the far north, where the red pixels dominate over the blue, you can see shadows of the features offset. E.g., Greenland is plotted twice, as is every other feature. Greenland is plotted were it is (fuzzy) but also over/north of Norway. So the bins computed from Step 1 are the wrong bins since the equatorial regions, mostly blue in the hits, are mostly from Step 2, and that maps looks mostly good. This is why with -n+a it looks good since that bypasses step 1. I will debug step 1 with tiles and with single grid and presumably I will learn something. I think I am getting closer.

[PaulWessel]

Above I am referring to features in the main plot, not the hitmap.

[PaulWessel]

Pretty sure it is because the assembly of the tiles passes in region as 0/360 but the plot is -60/300, and that is the ghost shift etc we are seeing. Jeez...

[PaulWessel]

Think I have a solution now. Just running tests but works for the example herein. Also think it points to a solution for the pyGMT case.

[PaulWessel]

It's 01d full-grid data, but is passed to GMT as an xarray.

Remind me (again), passes as GMT_VIA_MATRIX and REFERENCE or DUPLICATE?

[seisman]

• GMT_IN|GMT_IS_REFERENCE: Global earth relief data loaded as xarray doesn't work if projection center is not 0 pygmt#515 (comment)
• GMT_IN: Global earth relief data loaded as xarray doesn't work if projection center is not 0 pygmt#515 (comment)
• GMT_IN|GMT_IS_DUPLICATE: crash

[seisman]

FYI, #3799 doesn't fix GenericMappingTools/pygmt#515.

[PaulWessel]

No, it can't, but I know why. More tomorrow.

[seisman]

The issue was fixed in #3799. Closing.