[R] netCDF to raster and spatial projection
Frede Aakmann Tøgersen
frtog at vestas.com
Wed Apr 16 12:25:37 CEST 2014
Hi Bea
Well the first hit lead me to rasterProject{raster}. Will this suit you?
> rasterMG.proj <- projectRaster(rasterMG, crs=CRS("+init=epsg:21781"))
> print(rasterMG.proj)
class : RasterLayer
dimensions : 116, 91, 10556 (nrow, ncol, ncell)
resolution : 40.1, 40.1 (x, y)
extent : 478794.9, 482444, 646645.8, 651297.4 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:21781 +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs
data source : in memory
names : part.a
values : 0, 1 (min, max)
The difference can be seen by plotting.
plot(rasterMG)
plot(rasterMG.proj)
Yours sincerely / Med venlig hilsen
Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling
Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com
Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On Behalf Of Frede Aakmann Tøgersen
> Sent: 16. april 2014 12:11
> To: Bea GD; r-help at r-project.org
> Subject: Re: [R] netCDF to raster and spatial projection
>
> Well no need to have your data. Usually one can find some suitable data in
> the help for the functions under question. So here is a reproducible example
> from ?meuse.grid.
>
> > data(meuse.grid)
> > coordinates(meuse.grid) <- ~x+y
> > proj4string(meuse.grid) <- CRS("+init=epsg:28992") # see ?meuse for this
> crs
> > gridded(meuse.grid) <- TRUE
> > summary(meuse.grid)
>
> Object of class SpatialPixelsDataFrame
> Coordinates:
> min max
> x 178440 181560
> y 329600 333760
> Is projected: TRUE
> proj4string :
> [+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555
> +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
> +ellps=bessel +units=m +no_defs]
> Number of points: 3103
> Grid attributes:
> cellcentre.offset cellsize cells.dim
> x 178460 40 78
> y 329620 40 104
> Data attributes:
> part.a part.b dist soil ffreq
> Min. :0.0000 Min. :0.0000 Min. :0.0000 1:1665 1: 779
> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.1193 2:1084 2:1335
> Median :0.0000 Median :1.0000 Median :0.2715 3: 354 3: 989
> Mean :0.3986 Mean :0.6014 Mean :0.2971
> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.4402
> Max. :1.0000 Max. :1.0000 Max. :0.9926
>
> > rasterMG <- raster(meuse.grid)
> > print(rasterMG)
>
> class : RasterLayer
> dimensions : 104, 78, 8112 (nrow, ncol, ncell)
> resolution : 40, 40 (x, y)
> extent : 178440, 181560, 329600, 333760 (xmin, xmax, ymin, ymax)
> coord. ref. : +init=epsg:28992 +proj=sterea +lat_0=52.15616055555555
> +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
> +ellps=bessel +units=m +no_defs
> data source : in memory
> names : part.a
> values : 0, 1 (min, max)
>
> > rasterMG.proj <- spTransform(rasterMG, CRS("+init=epsg:21781"))
> Error in spTransform(rasterMG, CRS("+init=epsg:21781")) :
> load package rgdal for spTransform methods
> >
>
> Now perhaps doing the projection before casting it to a raster will work (yes
> I'm guessing since error thrown above is not very informative and traceback()
> does not give anything useful).
>
> > meuse.grid.proj <- spTransform(meuse.grid, CRS("+init=epsg:21781"))
> Warning message:
> In spTransform(meuse.grid, CRS("+init=epsg:21781")) :
> Grid warping not available, coercing to points
>
> > summary(meuse.grid.proj)
>
> Object of class SpatialPointsDataFrame
> Coordinates:
> min max
> x 479029.8 482197.1
> y 646927.4 651003.9
> Is projected: TRUE
> proj4string :
> [+init=epsg:21781 +proj=somerc +lat_0=46.95240555555556
> +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel
> +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs]
> Number of points: 3103
> Data attributes:
> part.a part.b dist soil ffreq
> Min. :0.0000 Min. :0.0000 Min. :0.0000 1:1665 1: 779
> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.1193 2:1084 2:1335
> Median :0.0000 Median :1.0000 Median :0.2715 3: 354 3: 989
> Mean :0.3986 Mean :0.6014 Mean :0.2971
> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.4402
> Max. :1.0000 Max. :1.0000 Max. :0.9926
>
> > gridded(meuse.grid.proj) <- TRUE
> suggested tolerance minimum: 0.737332
> Error in points2grid(points, tolerance, round) :
> dimension 1 : coordinate intervals are not constant
> >
>
> A warning and an error indicates that the projection results in something that
> is not on a regular grid.
>
> I don't know what to do but to read some more of the documentation for sp,
> rgdal, etc.
>
> Hopefully somebody comes up with something that helps us.
>
> Yours sincerely / Med venlig hilsen
>
>
> Frede Aakmann Tøgersen
> Specialist, M.Sc., Ph.D.
> Plant Performance & Modeling
>
> Technology & Service Solutions
> T +45 9730 5135
> M +45 2547 6050
> frtog at vestas.com<mailto:frtog at vestas.com>
> http://www.vestas.com<http://www.vestas.com/>
>
> Company reg. name: Vestas Wind Systems A/S
> This e-mail is subject to our e-mail disclaimer statement.
> Please refer to
> www.vestas.com/legal/notice<http://www.vestas.com/legal/notice>
> If you have received this e-mail in error please contact the sender.
>
> From: Bea GD [mailto:aguitatierra at hotmail.com]
> Sent: 16. april 2014 11:00
> To: Frede Aakmann Tøgersen; r-help at r-project.org
> Subject: Re: [R] netCDF to raster and spatial projection
>
> Hi Frede,
>
> Thanks so much for your reply!
>
> I've tried what you said but I get the following error:
>
>
>
> Error in spTransform(rasterDF1, crs("+init=epsg:21781")) :
>
> load package rgdal for spTransform methods
>
> I've checked search() and rgdal is before sp.
>
>
> > search()
>
> [1] ".GlobalEnv" "package:chron" "package:sm" "package:rgeos"
>
> [5] "package:maptools" "package:ncdf" "package:rgdal"
> "package:raster"
>
> [9] "package:sp" "tools:rstudio" "package:stats" "package:graphics"
>
> [13] "package:grDevices" "package:utils" "package:datasets"
> "package:methods"
>
> [17] "Autoloads" "package:base"
>
> Also, when I do library("rgdal") I get this message:
>
>
> > library("rgdal")
>
> rgdal: version: 0.8-16, (SVN revision 498)
>
> Geospatial Data Abstraction Library extensions to R successfully loaded
>
> Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
>
> Path to GDAL shared files: C:/Users/bgonzale/Documents/R/win-
> library/3.0/rgdal/gdal
>
> GDAL does not use iconv for recoding strings.
>
> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
>
> Path to PROJ.4 shared files: C:/Users/bgonzale/Documents/R/win-
> library/3.0/rgdal/proj
>
> Maybe the problem is something to do with this? I don't know how to fix it.
>
> I'd like to provide you with some data, how would be the best way to
> post/sahre a raster?
>
> Thanks!
>
>
> On 16.04.2014 10:17, Frede Aakmann Tøgersen wrote:
>
> Hi Beatriz
>
>
>
> Try to skip this step
>
>
>
> # Reprojecting into CH1903_LV03
>
> # First, change the coordinate reference system (crs)
>
> proj4string(rasterDF1) <- "+init=epsg:21781"
>
>
>
>
>
> And just do this
>
>
>
> # Second, reproject the raster
>
> rasterDF1.CH <- spTransform(rasterDF1, crs("+init=epsg:21781"))
>
>
>
>
>
> Also there is a spTransform both in the rgdal and sp packages. So are they
> masking each other? rgdal should be before sp in the search() list.
>
>
>
> I cannot be of more help since you provided no data.
>
>
>
>
>
> Yours sincerely / Med venlig hilsen
>
>
>
>
>
> Frede Aakmann Tøgersen
>
> Specialist, M.Sc., Ph.D.
>
> Plant Performance & Modeling
>
>
>
> Technology & Service Solutions
>
> T +45 9730 5135
>
> M +45 2547 6050
>
> frtog at vestas.com<mailto:frtog at vestas.com>
>
> http://www.vestas.com
>
>
>
> Company reg. name: Vestas Wind Systems A/S
>
> This e-mail is subject to our e-mail disclaimer statement.
>
> Please refer to
> www.vestas.com/legal/notice<http://www.vestas.com/legal/notice>
>
> If you have received this e-mail in error please contact the sender.
>
>
>
> -----Original Message-----
>
> From: r-help-bounces at r-project.org<mailto:r-help-bounces at r-project.org>
> [mailto:r-help-bounces at r-project.org]
>
> On Behalf Of Beatriz R. Gonzalez Dominguez
>
> Sent: 16. april 2014 08:22
>
> To: r-help at r-project.org<mailto:r-help at r-project.org>
>
> Subject: [R] netCDF to raster and spatial projection
>
>
>
> I've recently started using R for spatial data. I'd be really grateful
>
> if you could could help me with this. Thanks!
>
>
>
> Sorry I don't provide a reproducible example. Please ask me if you have
>
> any questions about the data.
>
>
>
> 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!
>
>
>
> ______________________________________________
>
> R-help at r-project.org<mailto:R-help at r-project.org> mailing list
>
> https://stat.ethz.ch/mailman/listinfo/r-help
>
> PLEASE do read the posting guide http://www.R-project.org/posting-
>
> guide.html
>
> and provide commented, minimal, self-contained, reproducible code.
>
>
>
>
>
>
> [[alternative HTML version deleted]]
More information about the R-help
mailing list