CODEWITHSUNDEEP
rmappingnetcdftemperature
mp_ylee
mp_ylee

Reputation: 41

Mapping temperature data from an .nc file in R

I downloaded temperature data from [NARR] (https://www.esrl.noaa.gov/psd/data/gridded/data.narr.monolevel.html) specifically "Air temperature at 2m" -monthly mean

I opened the file using the package "ncdf4". The data has 4 dimensions- time, x, y, nbnds. y corresponds to lat and x corresponds to lon. There is a variable (not dimension) called air which I do not know how to use, although this is the temperature information.

My end goal is to map the temperature data on a map of North America, using averaged temperature data for each month for each year (12 maps, one for each month).

I am having trouble identifying how to use the data as all of the dimensions are just really long lists of numbers that don't seem to have meaning (eg. the x coordinates look like this: 6232896 6265359 6297822 6330285 6362748 6395211 6427674 6460137 6492600 6525063 6557526 6589989, and so do the y values and time).

Here is the code I am using to view the dimensions:

temp2m <- nc_open("air.2m.mon.mean.nc")
time <- temp2m$dim$time$vals
lat <- temp2m$dim$x$vals
lon <- temp2m$dim$y$vals
nbnds <- temp2m$dim$nbnds$vals

If someone could help me view the data as well as map temperature data onto North America that would be great.

Thank you!

Upvotes: 2

Views: 1380

Answers (1)

Spacedman
Spacedman

Reputation: 94182

You can use the raster package to read these into a stack:

> library(raster)
> air = stack("./air.2m.mon.mean.nc")

(Note, you may need a raster package compiled with netcdf drivers...)

You can then plot them by slice or by time-name:

> plot(air[[23]])
> plot(air[["X1979.10.01.01.01.15"]])

enter image description here

The stack prints like this:

> air
class       : RasterStack 
dimensions  : 277, 349, 96673, 450  (nrow, ncol, ncell, nlayers)
resolution  : 32462.99, 32463  (x, y)
extent      : -16231.49, 11313351, -16231.5, 8976020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +x_0=5632642.22547 +y_0=4612545.65137 +lat_0=50 +lon_0=-107 +lat_1=50 +lat_2=50 +ellps=WGS84 
names       : X1979.01.01.00.01.15, X1979.02.01.00.01.15, X1979.03.01.00.01.15, X1979.04.01.01.01.15, X1979.05.01.01.01.15, X1979.06.01.01.01.15, X1979.07.01.01.01.15, X1979.08.01.01.01.15, X1979.09.01.01.01.15, X1979.10.01.01.01.15, X1979.11.01.00.01.15, X1979.12.01.00.01.15, X1980.01.01.00.01.15, X1980.02.01.00.01.15, X1980.03.01.00.01.15, ...

and those coordinates are not really lat-long, but are in a transformed coordinate system described by that "coord. ref." string. If you want to put it on a lat-long map you need to warp it:

> air_ll = projectRaster(air[[1]],crs="+init=epsg:4326")
> plot(air_ll)

enter image description here

It might be better for you to transform any other data to this system, and keep the grid unprojected. Just look up how to deal with spatial data in R for more info on projections and transformations.

Upvotes: 3

Related Questions