getting a set of topographic data

[mdsumner]

Here's Gebco 2014 to polar GeoTIFF

## key here is a_srs because the NetCDF doesn't have that, and gdalwarp doesn't complain!
gdal_translate -projwin -180 -20 180 -90 /rdsi/PRIVATE/raad/data_local/www.bodc.ac.uk/gebco/GEBCO_2014_2D.nc gebcosouth.vrt -of VRT -a_srs "+init=epsg:4326"

gdalwarp gebcosouth.vrt -te -8680000 -8680000 8680000 8680000 -tr 1000 1000 -t_srs "+proj=stere +lat_0=-90 +lat_ts=-71 +datum=WGS84 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0" gebco_ps.vrt -of VRT -r bilinear

gdal_translate gebco_ps.vrt gebco_ps.tif -co COMPRESS=LZW -co TILED=YES

Read and plot that fast with vapour

f <- "/mnt/temp/gebco/gebco_ps.tif"
ss <- vapour:::ssraster(f)
plot(vapour:::pull_ssraster(ss))

See latitudecircle and latmask in spex for more.

[mdsumner]

The gdalwarp approach doesn't actually work, you always get a seam at the dateline - Lucy struggled with this too. It's probably best to use projectRaster rather than gdalwarp, or some DIY trick.

[mdsumner]

The R approach.

The key is to not have a target pixel straddling the dateline, because it's centre pixel will fall exactly on the dateline. So, ensure the offset/resolutions combinations have an edge at target X = 0 (for lon_0=0 case).

library(raadtools)
tfile <- "C:/data/gebcosouth.grd"
if (file.exists(tfile)) {
  topo <- raster(tfile)
} else {
  topo <- readtopo("gebco_14")
  topo <- crop(topo, extent(-180, 180, -90, -20), filename = tfile)
}
topo <- readAll(topo)

resolutions <- c(16000, 4000, 1000)
filenames <- sprintf("southpolar_gebco14_%ikm.tif", resolutions/1000)
## make sure this is TRUE so we don't sample from the dateline exactly
offset <- 8688000
all(((offset / resolutions) - (offset %/% resolutions)) == 0)
prj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +datum=WGS84 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0"

for (i in seq_along(resolutions)) {
target1 <- raster(extent(-offset, 0, -offset, 0), res = resolutions[i], crs = prj)
target2 <- raster(extent(-offset, 0, 0, offset), res = resolutions[i], crs = prj)
target3 <- raster(extent(0, offset, 0, offset), res = resolutions[i], crs = prj)
target4 <- raster(extent(0, offset, -offset, 0), res = resolutions[i], crs = prj)
r1 <- projectRaster(topo, target1)
r2 <- projectRaster(topo, target2)
r3 <- projectRaster(topo, target3)
r4 <- projectRaster(topo, target4)
r <- merge(r1, r2, r3, r4)
writeRaster(r, filenames[i], options = c("COMPRESS=DEFLATE", "TILED=YES"), datatype = "INT2S")
}

[mdsumner]

The GDAL approach

try -wo options specifically SOURCE_EXTRA

https://gis.stackexchange.com/questions/198954/gdalwarp-reprojects-only-a-portion-of-the-raster

https://gis.stackexchange.com/questions/213698/gdalwarp-malfunctioning-for-geos-projection

http://www.gdal.org/gdalwarp.html

http://www.gdal.org/structGDALWarpOptions.html#a0ed77f9917bb96c7a9aabd73d4d06e08

[mdsumner]

SOmap covers this now