2
votes

I used the paleoview software to download some variables from past climate, including the mean temperature from 15000BP-10000BP (I could upload a file on request but its a GB at least).

The main problem is that when I read the raster, it contains only positive latitude and positive longitude. I know from the article that this has a 2.5*2.5 resolution.

Using Raster

I loaded both the raster and ncdf4 libraries to read it using raster

library(raster)
library(ncdf4)

When I read it using the following code

r <- raster("mean_temperature-15000BP-10000BP.nc",  varname = "14000BP-13000BP/13300BP")

I get the following info

r
class       : RasterLayer 
band        : 1  (of  12  bands)
dimensions  : 72, 144, 10368  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0.5, 144.5, 0.5, 72.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : Mean.Temperature 
z-value     : 1 
zvar        : 14000BP-13000BP/13300BP 

As you see the extent is only positive, but I know it has the data for the whole world, when I plot the map I get the following image:

plot(r)

enter image description here

which clearly shows the patterns expected with the polar circles having extremely low temperatures and the Antarctica being larger that the arctic.

clearly the +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 is wrong here, and I think that if I figure out whats the actual coord. ref. I could reproject it and get the raster in the right format

Using ncdf4

I tried working with the ncdf4 package to try to get more info from the layer, and this is what I did:

nc <- nc_open("mean_temperature-15000BP-10000BP.nc")

After reading the layer using the nc_open function I see the names of the variables to try to understand more of the layers, here I show the first 10

names(nc$var)[1:10]

[1] "window"                  "width"                   "decimals"                "months"                  "latitudes"               "longitudes"              "15100BP-15000BP/15100BP" "15100BP-15000BP/15099BP"
[9] "15100BP-15000BP/15098BP" "15100BP-15000BP/15097BP"

So if I keep looking and go to the latitude and longitude names I get:

ncatt_get(nc, attributes(nc$var)$names[5])
$units
[1] "degrees north"

and

ncatt_get(nc, attributes(nc$var)$names[6])
$units
[1] "degrees east"

Any idea on how to reproject this rasters to get a latitude that goes from -90, 90 and longitude from -180 to 180 that you would expect in a +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 projections

1

1 Answers

1
votes

To fix this, you could set the extent:

library(raster)
r <- raster("mean_temperature-15000BP-10000BP.nc",  varname = "14000BP-13000BP/13300BP")
extent(r) <- c(-180,180,-90,90)

But note that sometimes longitude in climate data goes from 0 to 360 (not here, I think):

extent(r) <- c(0,360,-90,90)
r <- rotate(r)

To see more about what's in the ncdf, you can do

print(r)

Of course you only want to do this is you are really sure about the extent, and you may want to contact the data provider to report this issue.