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." "????? ?????? ???? ?"
Extract spatial mean of subset of netCDF.
3 messages · Aseem Sharma, Ben, Robert J. Hijmans
13 days later
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
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