[R-sig-Geo] Extract spatial mean of subset of netCDF.
Robert J. Hijmans
r.hijmans at gmail.com
Sat Dec 20 07:51:07 CET 2014
Aseem,
The problem is that you are trying to subset a raster using
latitude/longitude coordinates, while the raster has another CRS.
You can fix this, I think, by assigning the correct crs (you can find
these at http://www.spatialreference.org/)
and then, for example, adjusting your extent, or by using projectRaster
Robert
library(raster)
nc1 <- "GlobSnow_SWE_L3A_19900101_v2.0.nc"
nc2 <- "GlobSnow_SWE_L3A_19900102_v2.0.nc"
nc_data <-stack(nc1, nc2, varname="SWE")
nc_data
crs(nc_data) <- "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0
+a=6371228 +b=6371228 +units=m +no_defs"
egeo <- extent(-78,-74,50,52)
x <- raster(egeo, crs='+proj=longlat +datum=WGS84' )
e <- extent(projectExtent(x, "+proj=laea +lat_0=90 +lon_0=0 +x_0=0
+y_0=0 +a=6371228 +b=6371228 +units=m +no_defs"))
d <- crop(nc_data, e)
cellStats(d, mean)
# or
as.vector(cellStats(d, mean))
### or
z <- projectRaster(nc_data, x)
cellStats(d, mean)
# number are similar, but (as expected) not the same.
Please provide more complete information, such as warning or error
messages, when you encounter a problem. Show the output of your script
(as below) that I needed to see to understand what is going on.
> library(raster)
> nc1 <- "GlobSnow_SWE_L3A_19900101_v2.0.nc"
> nc2 <- "GlobSnow_SWE_L3A_19900102_v2.0.nc"
> nc_data <-stack(nc1, nc2, varname="SWE")
Warning messages:
1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
spatial_ref=PROJCS["NSIDC EASE-Grid North",GEOGCS["Unspecified datum
based upon the International 1924 Authalic
Sphere",DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",SPHEROID["International
1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],AUTHORITY["EPSG","6053"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4053"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3408"]];
GeoTransform=-9036842.762 25067.525 0 9036842.763000002 0 -25067.525
2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
grid_mapping_name; false_easting; false_northing;
scale_factor_at_projection_origin; scale_factor_at_central_meridian;
standard_parallel; standard_parallel1; standard_parallel2;
longitude_of_central_meridian; longitude_of_projection_origin;
latitude_of_projection_origin; straight_vertical_longitude_from_pole;
longitude_of_prime_meridian; semi_major_axis; inverse_flattening;
+proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2; +lon_0; +lon_0;
+lat_0; +lon_0; +lon_0; +a; +rf
3: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
spatial_ref=PROJCS["NSIDC EASE-Grid North",GEOGCS["Unspecified datum
based upon the International 1924 Authalic
Sphere",DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",SPHEROID["International
1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],AUTHORITY["EPSG","6053"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4053"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3408"]];
GeoTransform=-9036842.762 25067.525 0 9036842.763000002 0 -25067.525
4: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
grid_mapping_name; false_easting; false_northing;
scale_factor_at_projection_origin; scale_factor_at_central_meridian;
standard_parallel; standard_parallel1; standard_parallel2;
longitude_of_central_meridian; longitude_of_projection_origin;
latitude_of_projection_origin; straight_vertical_longitude_from_pole;
longitude_of_prime_meridian; semi_major_axis; inverse_flattening;
+proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2; +lon_0; +lon_0;
+lat_0; +lon_0; +lon_0; +a; +rf
> nc_data
class : RasterStack
dimensions : 721, 721, 519841, 2 (nrow, ncol, ncell, nlayers)
resolution : 25067.52, 25067.52 (x, y)
extent : -9036843, 9036843, -9036843, 9036843 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
names : Daily.Snow.Water.Equivalent.1, Daily.Snow.Water.Equivalent.2
On Fri, Dec 19, 2014 at 8:55 AM, Benjamin.Galuardi <galuardi at mit.edu> wrote:
> Hi Aseem,
>
> I can offer the following advice:
>
> 1) there is no
>
> dt.ras.c
>
> variable created in your code. It would be helpful to include that line in
> the question.
>
> 2) do you need a raster brick? unless you have very large data, this step
> seems unnecessary.
> 3) are you looking for a mean across layers or a mean for each layer within
> the stack?
> 4) try using 'extract'
> mean( extract(raster, extent), na.rm=T)
> 4a) use na.rm=T! I'm not sure but perhaps this is causing problems.
> 5) double check the import. if the orientation of the raster is not what you
> think it is, you ,ay be extracting the wrong area.
>
> hope that helps
>
> Ben
>
>
> On 12/6/2014 11:47 AM, Aseem Sharma wrote:
>>
>> Hi,
>>
>> I am new with raster work in R.
>>
>> I am trying to get the spatial mean of a subset (define by the extent) of
>> larger domain netCDF data file using raster and ncdf packages
>>
>> The netCDF data of Snow Water Equivalent(SWE) were obtained from
>> http://www.globsnow.info/swe/archive_v2.0/. The two files I have used are
>> nc1
>>
>> <https://www.dropbox.com/s/tk388w67tjg704s/GlobSnow_SWE_L3A_19900101_v2.0.nc?dl=0>
>> and nc2
>>
>> <https://www.dropbox.com/s/pft2dlio96n6ncw/GlobSnow_SWE_L3A_19900102_v2.0.nc?dl=0>
>> .
>>
>> The codes I used are
>>
>> library(raster)
>>
>> library(ncdf)
>>
>> nc1<-"C:\\Users\\bunug_000.BUNU\\Documents\\Dropbox\\Public\\
>> GlobSnow_SWE_L3A_19900101_v2.0.nc"
>>
>> nc2<-"C:\\Users\\bunug_000.BUNU\\Documents\\Dropbox\\Public\\
>> GlobSnow_SWE_L3A_19900102_v2.0.nc"
>>
>>
>>
>> about_nc<-open.ncdf(nc1)
>>
>> print(about_nc)
>>
>>
>>
>> nc_data<-stack(nc1,nc2,varname="SWE")
>>
>> print(nc_data)
>>
>> dim(nc_data)
>>
>> nc_data1<-brick(nc_data)
>>
>> plot(nc_data1)
>>
>> subsetd<-extent(-78,-74,50,52)
>>
>> subsetd
>>
>> crop_data <- crop(nc_data1, subsetd)
>>
>> dim(dt.ras.c)
>>
>>
>>
>> dt.anm<-as.data.frame(cbind(unname(cellStats(dt.ras.c,
>> stat="mean")),seq(1:2))) # use unname to separate named numebr to just
>> numbers
>>
>> str(dt.anm)
>>
>> colnames(dt.anm)<-c("Mean_SWE","ID")
>>
>>
>> What is not right here so that I am not getting mean SWE. It should not be
>> -1.
>>
>> Any suggestions will be highly appreciated.
>>
>> Thank you.
>>
>>
>> ------------------
>>
>> "Namaste नमस्ते"
>>
>> Aseem Sharma
>>
>> Graduate Research Assistant
>>
>> Northern Hydrometeorology Group(NHG)
>>
>> Natural Resources and Environmental Studies Institute(NRESi)
>>
>> University of Northern British Columbia
>>
>> Prince George, BC, V2N 4Z9, Canada
>>
>> Tel: 250-960-5427
>>
>> Web: http://www.unbc.ca/
>>
>>
>> "All know the Way, but few actually walk it."
>> "सबैको कल्याण होस् ।"
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> --
> ==================================================
> Benjamin Galuardi
> PhD Student
> NMFS-Seagrant Population Dynamics Fellow
> SMAST, UMass Dartmouth
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
More information about the R-sig-Geo
mailing list