Question

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!

Was it helpful?

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"))

OTHER TIPS

Here is another option GDAL. Using gdal_translate you can convert netcdf file with bands to geotiffs.

 gdal_translate -a_srs EPSG:4326 NETCDF:File_Name.nc:Band_Name  -of ‘Gtiff’ Output_FileName.geotiff

To explore more options in gdal_translate you can visit this link

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top