[gdal-dev] downsampling geotiff with a low-pass filter

[JDA]

> It's better to use a sensible smoothing method in the first place. I 
> suppose noone has given this over much thought as in the past you were 
> ever so happy about every bit of resolution you could get. But in a time 
> where we get very high resolution LIDAR data, the need to downsample 
> properly is arising. Look at the interpolation methods: gdalwarp lists 
> twelve different ones. The first few are for upsampling, and the 
> remainder mostly for dealing with noisy data. Upsampling is well 
> covered: cubicspline and lanczos are reasonably sophisticated upsampling 
> filters, but there is no good downsampling filter. I think this is an 
> omission, hence my post. The problem is real; downsampling with 
> 'average' produces artifacts, even from previously upsampled data.

There’s a wide array of smoothing options available If you’re willing to work in python. Based on https://gis.stackexchange.com/a/10467, the basic idea is to load the raster into a numpy array and then convolve it with either a kernel of arbitrary size. 

I’ve written a program implementing that idea with both a Gaussian and a mean kernel here:
https://github.com/jacobdadams/general_scripts/blob/master/raster_chunk_processing.py. It only smooths the data without downsampling, but there may be some python functions that downsample as well. 

There’s a lot of code in there for dealing with massive rasters by processing them in chunks in parallel, but the blur_mean and blur_gauss functions are where the smoothing is done. I’ve written an installation and usage guide at https://gisjake.blogspot.com/2018/10/rasterchunkprocessingpy-installation.html?m=1. 

I’ve successfully used both the Gaussian blur and the Mesh Denoise methods to get rid of the surface “noise” (grass, alfalfa fields, etc) on half-meter lidar. I’ve combined the resulting hillshade with a 10 meter hillshade using gdalbuildvrt, but I haven’t found a good way to combine the 10m and half-meter DEMs yet.  

The valley floor slope layers on the Cache County parcel and zoning viewer were produced from the smoothed lidar DEM (cachecounty.org/gis, click on Parcel and Zoning Viewer, or direct link: http://66.232.67.238/Websites/Parcel%20and%20Zoning%20Viewer/; in the map, open the layer list, expand Sensitive/Planning Layers, turn on the two slope layers, then zoom in close so the parcel numbers don’t cover everything). 

Jake
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20181116/8ff0511b/attachment.html>