[R] Problems trying to place a global map with Ncdf data plot
Bert Gunter
bgunter@4567 @end|ng |rom gm@||@com
Sun Feb 17 16:38:18 CET 2019
The r-sig-geo list would probably be a better place to post this, as they
specialize in this sort of thing.
Bert Gunter
"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
On Sun, Feb 17, 2019 at 3:48 AM rain1290--- via R-help <r-help using r-project.org>
wrote:
> Hello there,
>
> I am trying to overlay a global map with ncdf data of precipitation for a
> particular location (using specific coordinates). The file is in ncdf
> format
> (commonly used to store away climate data), and I am currently attempting
> to
> place a global map on plotted precipitation values. However, I am having
> difficulty placing a global map on this plot and am encountering errors. I
> will show you what I have done:
>
> #To create a plot of precipitation data using the following ncdf file - the
> following works fine and provides the distributions global precipitation
> values (Land+Water values):
>
> library(ncdf4)
> Can<-"MaxPrecCCCMACanESM2rcp45.nc"
>
>
> >Model<-nc_open(Can)
> >print(Model)
> >attributes(Model$var)
> >$names
> >dat<-ncvar_get(Model, "onedaymax")
> >dat[128,50,1] #View onedaymax for selected latitude, longitude and Year
> >nc_lat<-ncvar_get(Model,attributes(Model$dim)$names[2]) #Retrieve latitude
> >nc_lon<-ncvar_get(Model,attributes(Model$dim)$names[3]) #Retrieve
> longitude
> >print(paste(dim(nc_lat), "latitudes and", dim(nc_lon), "longitudes"))
> >library(maptools)
> >map<-dat[,,5] #Precipitation for all longitudes, latitudes, and Year 5
> >grid<-expand.grid(nc_lon=nc_lon, nc_lat=nc_lat)
> >image(nc_lon,nc_lat,map, ylab="Latitude", xlab="Longitude", main="One-day
> maximum precipitation")
> >levelplot(map~nc_lon*nc_lat, data=grid, at=cutpoints, cuts=11,
> ylab="Latitude", xlab="Longitude", >main="Year 5 one-day maximum
> precipitation (mm/day) for CanESM2 under RCP4.5", pretty=T,
> col.regions=(rev(brewer.pal(10, "Spectral"))))
>
> #To place a global map on the map that map that returns using the above
> code. *This is where errors begin:
>
> >ggplot()+geom_point(aes(x=nc_lon,y=nc_lat,color="onedaymax"),
> size=0.8)+borders("world",
>
> colour="black")+scale_color_viridis(name="onedaymax")+theme_void()+coord_quickmap()
> *Error: Aesthetics must be either length 1 or the same as the data (128):
> x,
> y, colour*
>
>
> Why doesn't this work? Could it be that I am not including the "time"
> dimension in the ggplot function? The problem, though, is when I try to
> obtain the "time" dimension, like I did for latitude and longitude, I
> receive the following error:
>
> t<-ncvar_get(Model,"time")
> *Error in nc$dim[[idobj$list_index]] :
> attempt to select more than one element*
>
> If it helps, this is what the variables and dimensions look like in the
> ncdf
> file:
>
> /File MaxPrecCCCMACanESM2rcp45.nc (NC_FORMAT_NETCDF4):
>
> 3 variables (excluding dimension variables):
> double onedaymax[lon,lat,time] (Contiguous storage)
> units: mm/day
> double fivedaymax[lon,lat,time] (Contiguous storage)
> units: mm/day
> short Year[time] (Contiguous storage)
>
> 3 dimensions:
> time Size:95
> lat Size:64
> units: degree North
> lon Size:128
> units: degree East/
>
> Any help would be greatly appreciated!!!!
>
> Thanks,
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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