Search code examples
projectionrastercoordinate-systemsnetcdf

netCDF to raster and projection


I've recently started using R for spatial data. I'd be really grateful if you could could help me with this. Thanks!

I've extracted data from a multidimensional netCDF file. This file had longitude, latitude and temperature data (for 12 months of a specific year).

From this netCDF I've got a data frame for January with these variables: longitude, latitude, temperature.

With this data frame I've created a raster.

# Packages
library("sp")
library("raster") 
library("rgdal")
library("ncdf")
library("maptools")
library("rgeos")
library("sm")
library("chron")

# Dataframe to raster
# Create spatial points data frame
coordinates(tmp.df01) <- ~ lon + lat
# Coerce to SpatialPixelsDataFrame
gridded(tmp.df01) <- T
# Coerce to raster
rasterDF1 <- raster(tmp.df01)
> print(tmp.df01)
class       : SpatialPixelsDataFrame 
dimensions  : 103, 241, 24823, 1  (nrow, ncol, npixels, nlayers)
resolution  : 0.02083333, 0.02083333  (x, y)
extent      : 5.739583, 10.76042, 45.73958, 47.88542  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
names       :           TabsM_1 
min values  : -18.1389980316162 
max values  :  2.26920962333679 

There is no value for 'coord. ref.'

The projection of the original netCDF was WGS84. So I gave this projection to the raster.

proj4string(rasterDF1) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Then, I wanted to reproject my raster to another projection:

# Reprojecting into CH1903_LV03
# First, change the coordinate reference system (crs)
proj4string(rasterDF1) <- "+init=epsg:21781"
# Second, reproject the raster
rasterDF1.CH <- spTransform(rasterDF1, crs("+init=epsg:21781"))

At this point I get the following error:

Error in spTransform(rasterDF1, crs("+init=epsg:21781")) : 
  load package rgdal for spTransform methods

But the package rgdal is already uploaded! It must be something wrong in the code!


Solution

  • Here the code to solve the problem described.

    Solution provided by Frede Aakmann Tøgersen.

    tmp.df01 # tmp.df01 is a data.frame
    coordinates(tmp.df01) <- ~ lon + lat # tmp.df01 is now a SpatialPointsDataFrame
    
    # Assign orignial data projection
    proj4string(tmp.df01) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
    
    gridded(tmp.df01) <- T # tmp.df01 is now a SpatialPixelFrame
    
    # Coerce to raster
    rasterDF1 <- raster(tmp.df01) # rasterDF1 is a RasterLayer
    
    # To reproject the raster layer
    rasterDF1.proj <- projectRaster(rasterDF1, crs=CRS("+init=epsg:21781"))