Task #552
openFacilitate raster buffering across antimeridian
0%
Description
For DEM smoothing, along with other procedures used in the layer production pipeline, we need one or more of the raster inputs to have an extent that is slightly larger than the extent of the desired target region itself. Even if not strictly necessary for the routine to run without error, this spatial buffering helps avoid creation of edge artifacts resulting from running routines on individual tiles (i.e., separate spatial subsets), as we frequently do for computational or other reasons.
Although it is usually straightforward to ensure that the input rasters are sufficiently large, this is more complicated at the antimeridian, which serves as both the left (180W lon) and right (180E lon) extent of some key global raster inputs, and which is consequently associated with an artificial break in the real-world contiguity across this line. Relatively little terrestrial area spans this meridian, but there is some. It would be convenient if our tools were all aware that e.g. a request for column of pixels beyond 180E (as would be needed to buffer tiles along the eastern edge of the global extent) is in fact simply a request for the column of pixels adjacent 180W. But given that this isn't the case, we could benefit from having a consistent general approach for dealing with this problem independent of the particular procedure.
Note that this is an issue regardless of whether we're operating on the entire global layer at once or operating on tiles adjacent/crossing the antimeridian.
Updated by Jim Regetz over 11 years ago
One approach, hopefully useful at least in some cases, would be to create VRTs that (virtually) represent raster data repeated beyond the global extent in both the east and west directions, by whatever amount is deemed necessary. These VRTs would reference their locations in the way we need -- i.e., using lat coordinates <-180 degrees or >180 degrees -- but for the raster values themselves would simply point to the corresponding window within the strict global extent as stored in the original raster file(s).
For additional convenience, we could then create another single "global-plus" VRT that provides a single coherent view of the entire global dataset plus the east-west extensions. Any given processing routine could then use GDAL to subselect from this global+ VRT to generate whatever input is desired.
Illustration with a couple of GDEM2 tiles located on atlas:
# point to directory of resampled GDEM2 tiles for demo purposes
GDEM2DIR=/home/layers/data/terrain/dem-aster-gdem2-90m-resamp
# note that these two tiles are on opposite ends of the global extent,
# but contiguous in the real world:
gdalinfo $GDEM2DIR/N65E179_ReSample.tif
gdalinfo $GDEM2DIR/N65W180_ReSample.tif
# create new file that expresses the eastern tile as an extension of the
# western edge, effectively repeating this tile but cheaply as a VRT
gdal_translate -of VRT -a_ullr -181 66 -180 65 \
$GDEM2DIR/N65E179_ReSample.tif N65W181_ReSample.vrt
# now build VRT that captures all of these
gdalbuildvrt gdem2_extended.vrt \
$GDEM2DIR/N65{W180,E179}_ReSample.tif N65W181_ReSample.vrt
# generate sample 1-degree-wide tif tile centered on the antimeridian
gdal_translate -projwin -180.5 66 -179.5 65 gdem2_extended.vrt \
aster_test.tif