projectRaster: Raster projection of bathymetry data ("NOAA.nc') in the pacific – worldmap with marmap in the pacific

Question

I would like to reproduce the last example from the marmap vignette: 'marmap-DataAnalysis' for a pacific region. The example shows the orthographic projection of the world centered at lon = 50. Here is the example:

library(marmap)
library(raster)
# Get data for the whole world. Careful: ca. 21 Mo!
world <- getNOAA.bathy(-180, 180, -90, 90, res = 15, keep = TRUE)

# Switch to raster
world.ras <- marmap::as.raster(world)

# Set the projection and project
my.proj <-   "+proj=ortho +lat_0=0 +lon_0=50 +x_0=0 +y_0=0"
world.ras.proj <- projectRaster(world.ras,crs = my.proj)

# Switch back to a bathy object
world.proj <- as.bathy(world.ras.proj)

# Set colors for oceans and land masses
blues <- c("lightsteelblue4", "lightsteelblue3",
           "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# And plot!
plot(world.proj, image = TRUE, land = TRUE, lwd = 0.05,
     bpal = list(c(0, max(world.proj, na.rm = T), greys),
                 c(min(world.proj, na.rm = T), 0, blues)),
     axes = FALSE, xlab = "", ylab = "")
plot(world.proj, n = 1, lwd = 0.4, add = TRUE)

However, I would like to change the central to a pacific meridian, e.g. lon = 155.5. I tried this by changing the projection parameters to,

my.proj <- "+proj=ortho +lat_0=20 +lon_0=155.5 +x_0=0 +y_0=0"

but then,

world.ras.proj <- projectRaster(world.ras,crs = my.proj)

results in:

Error in if (nr != x@nrows | nc != x@ncols) { : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1],     xy[,  :
  259 projected point(s) not finite
2: In rgdal::rawTransform(projection(raster), crs, nrow(xy), xy[,     1],  :
  4 projected point(s) not finite

How can could I plot the 'bathymetric world' in a pacific region?

Robert Hijmans · Answer 1 · 2017-11-02T02:19:15.107

2

I have simplified your question (always good, and for me the data download did not work). In essence:

library(raster); library(rgdal)
prj1 <- "+proj=ortho +lat_0=0 +lon_0=0 +x_0=0 +y_0=0"
prj2 <- "+proj=ortho +lat_0=20 +lon_0=155.5 +x_0=0 +y_0=0" 
r <- raster()
r <- init(r, 'col')

# works
x1 <- projectRaster(r, crs = prj1)
# fails
x2 <- projectRaster(r, crs = prj2)

This is a bug. I have fixed it in raster version 2.6-2 (under development, should be available next week or so)

edited Nov 02 '17 at 02:19

answered Nov 02 '17 at 01:37

Robert Hijmans

40,301
4
55
63

Thank you. Looking forward to the update. Any more details on this bug? Is it related to the type of projection (+proj=ortho)? – Tony Nov 03 '17 at 00:14
The problem was with projection across the date-line – Robert Hijmans Nov 21 '17 at 18:45

Benoit · Answer 2 · 2017-11-21T20:52:35.737

This can be solved in marmap with the current/previous version of the raster package. You have to use the antimeridian=TRUE argument of the getNOAA.bathy() function and some trickery to allow the computations of the projection by the raster package.

The first trick is to download data with lon1 = lon2 = 0 since the antimeridian downloads 2 distinct datasets: from the antimeridian to lon1 and from lon2 to the antimeridian. Setting lon1 and lon2 equal to 0 downloads the whole world.

Then, you have to manually switch back to values of longitudes between -180 and 180 (and not 0 to 360 as produced by the animeridian argument of getNOAA.bathy()), hence the rownames(world2) <- ... line.

Finally, you have to apply the same -180 correction to specify the projection. Here is the code:

library(marmap)
library(raster)
# Get data for the whole world. Careful: ca. 21 Mo!
world2 <- getNOAA.bathy(0, 0, -90, 90, res = 15, keep = TRUE, antimeridian=TRUE)
rownames(world2) <- as.numeric(rownames(world2))-180

# Switch to raster
world.ras <- marmap::as.raster(world2)

# Set the projection and project
my.proj <-   "+proj=ortho +lat_0=20 +lon_0=155-180 +x_0=0 +y_0=0"
world.ras.proj <- projectRaster(world.ras,crs = my.proj)

# Switch back to a bathy object
world.proj <- as.bathy(world.ras.proj)

# Set colors for oceans and land masses
blues <- c("lightsteelblue4", "lightsteelblue3",
           "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# And plot!
plot(world.proj, image = TRUE, land = TRUE, lwd = 0.05,
     bpal = list(c(0, max(world.proj, na.rm = T), greys),
                 c(min(world.proj, na.rm = T), 0, blues)),
     axes = FALSE, xlab = "", ylab = "")
plot(world.proj, n = 1, lwd = 0.4, add = TRUE)

And here is the result:

The new version of raster (2.6-7) solves the problem of projecting across the date-line. However, due to rounding errors when downloading bathymetric data from NOAA servers, some missing cells may appear in plots. Here is an example with the code you posted in your original question:

And here is the summary() of the data:

summary(world)
# Bathymetric data of class 'bathy', with 1440 rows and 720 columns
# Latitudinal range: -89.88 to 89.88 (89.88 S to 89.88 N)
# Longitudinal range: -179.88 to 179.88 (179.88 W to 179.88 E)
# Cell size: 15 minute(s)

# Depth statistics:
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  -10635   -4286   -2455   -1892     214    6798 
# 
# First 5 columns and rows of the bathymetric matrix:
#          -89.875 -89.625 -89.375 -89.125 -88.875
# -179.875    2746    2836    2893    2959    3016
# -179.625    2746    2835    2892    2958    3015
# -179.375    2746    2835    2891    2957    3014
# -179.125    2746    2834    2890    2956    3013
# -178.875    2746    2834    2889    2955    3012

Hence, the solution using antimeridian=TRUE detailed above should be best.

projectRaster: Raster projection of bathymetry data ("NOAA.nc') in the pacific – worldmap with marmap in the pacific

2 Answers2