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)
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
print(nc)
too? – Tung