Downsampling SRTM15+V2.grd to other resolutions

[PaulWessel]

Our beta earth_relief_xxy grids were downsampled by spherical Gaussian filtering. This means that the wavelengths of features anywhere on earth where filtered the same way. However, the way the original grid was produced via Cartesian gridding of strips, there is much higher spectral resolution in the east-west direction as we go to higher latitudes, since the bins per degree longitude is fixed. Potentially, there is information at a shorter wavelength in east-west at high latitudes that is shorter than the pixel spacing at Equator.

My question then is this: Do we continue to do the spherical filtering so all parts of the planet has the same spectral content, or do we use a filter that does not reach wider in east-west as we go to high latitude? Below is a series of plots. The first three are from Africa near (0,0). The first plot is the source (15s) while the others are 5m smoothed versions. The second plot is spherically filtered while the third is Cartesian filtered. There are some tiny differences between the last two:

The next three are from NE Greenland at latitudes 75-80. Same order of plots and now you can see a clear difference in that the Cartesian filtering leaves more details in the east-west direction:

I don't really like that the original grid is not spatially uniform, but I cannot do much about that. I tend to like the spherical treatment I gave these grids, but you may feel otherwise. What do you guys think, @GenericMappingTools/core?

[PaulWessel]

To clarify how the spherical filtering works, see figures below. I show plots of the 15 arc second data distribution (small circles) and the desired output locations for the 2 minute grid (large black circles). The full filter size is plotted in pink (centered on (0,0) output node and has a full width 7.4 km equaling twice the desired grid spacing at Equator. This circle is the locus of all points on the grid that is within 7.2 km of the output node. Large circles indicate the filter areas for other output nodes. The first plot is for nodes near Equator:

Second plot is for nodes along 60N where the 7.4-km diameter circles cover more longitudinal range:

[PaulWessel]

Trying to get Walter to comment on this but he has yet to sign up for a GitHub account... Anyone else with comments? I think it would be nice to release v2 of these grids with GMT 6.

[seisman]

I don't know the technical details, but the spherical filtering makes more sense to me.

[PaulWessel]

Yes, me too, but I feel bad obliterating actual real data just because the grid has variable spatial resolution....

[joa-quim]

I think that the Cartesian makes more sense since the data originally has more resolution in the E-W direction. At sea, the situation may be even worst. Given the poorer bathymetry data coverage there are several places where nodes were filled by pure interpolation (I mean, no data control inside the cell). The spherical filtering will give more weight to those purely interpolated nodes (because it will reach more nodes in E-W).

[PaulWessel]

I will need to experiment but one problem by going Cartesian is the lack of proper boundary conditions at the pole and possibly at the periodic longitude borders. I need to see if we can combine the Cartesian radius with full geographic BCs.

[WalterHFSmith]

The problem, as I see it, boils down to whether you want to simply resample, or whether you really want to filter. Furthermore, by doing a "spherical" filter you are assuming that the filter impulse response ought to be isotropic, that is to say, independent of azimuth.
Seems to me that if one is going to treat latitude and longitude as Cartesian, meaning make a grid with sampling intervals dx and dy that are a fixed increment of lon and lat, then this grid already has a crazy sampling scheme. Suppose it is 15 arc-seconds, and you want to go to 1 arc-min with the same crazy sampling scheme. Then I think you want a block mean or block median of 4x4 blocks. OTOH, if you were trying to fit spherical harmonics to either of these (the 15 sec or the 1 min) then there would have to be some least-squares fitting because the spherical harmonics, even if carried to maximum degree and order 360x60 or 360x60x4, would necessarily smooth the high latitudes.
So it is a question of sampling versus smoothing in my mind.
Perhaps a compromise is to treat rectangularly sampled data as being resampled by a Fourier (sinc) reconstruction, followed by a low-pass rectangular filtering down to a wider sinc function.
If longitude is rectangularly sampled and spans 360 degrees it could be given a periodic boundary condition this way.
If latitude is sampled from -90 to 90 this way it is also possible to give it a periodic extension or appropriate boundary condition by considering it to extend across the poles and down the meridian on the other side of the Earth. I have had some luck doing pseudo-spherical filtering of global rectangular-sampled grids this way.

Even the sinc interpolant resampling has huge losses. See Fig 11 attached to my comment on pixel versus grid registrations.

[WalterHFSmith]

I should have added, the mis-use (in my opinion) of the word "resolution" is a pet peeve of mine. I think it would be helpful, in discussions like this, if we were to try to distinguish between "sampling" (or "sampling interval" or "sampling frequency", as necessary), which refers to the distance between samples, and not "resolution", which (in my view) properly refers to the scale of feature that stand out above the background noise. A rectangular sampling of lat and lon with N samples in latitude and 2N samples in longitude does not really "resolve" spherical harmonics to degree and order N.

[PaulWessel]

Yes I agree. Perhaps we can clean up some language in the use of that term as it relates to these files in 6.1.

[Esteban82]

I propose to send this for discussion. If we close it, I think that what is said here (which I think is important) would be hidden among the other closed issues.

[PaulWessel]

Sure, we can do that. I might add some more comments since in my mind "sampling" is to pick values at (presumably) a regular interval, but we actually do filtering first so down-sampling is always a mix of filtering and sampling the filtered grid, otherwise we would get nasty aliasing.