[R] netCDF to raster and spatial projection
Bea GD
aguitatierra at hotmail.com
Wed Apr 16 13:33:11 CEST 2014
Hi Frede,
Thanks so much! It seems to work perfectly, exactly what I wanted to do! :)
All the best,
Bea
On 16.04.2014 12:25, Frede Aakmann Tøgersen wrote:
> 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