[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