update

R/namer.R

@@ -1,5 +1,6 @@
 library(namer)
 
-filename <- "/Users/bartlein/Dropbox/DataVis/REarthSysSci/hdf5_intro.Rmd" # path to and name of .Rmd file 
+filename <- "/Users/bartlein/Dropbox/DataVis/REarthSysSci/plots1.Rmd" # path to and name of .Rmd file 
+unname_chunks(filename)
 name_chunks(filename)
 

R/plots_01.R

@@ -9,7 +9,6 @@ library(ggplot2)
 attach(sumcr)
 plot(Length)
 
-qplot(seq(1:length(Length)), Length)
 ggplot(data=sumcr, aes(x = Length)) + 
   geom_point(aes(x = 1:nrow(sumcr), y = Length)) + 
   geom_line(aes(x = 1:nrow(sumcr), y = Length)) +
@@ -69,6 +68,7 @@ ggplot(specmap, aes(x = O18)) +
   geom_line(stat = "density") +
   geom_rug(data = specmap, aes(x = O18)) 
 
+
 # scatter plots
 # use Oregon climate-station data 
 
@@ -147,7 +147,11 @@ detach(sumcr)
 plot(orstationc$lon, orstationc$lat, type="n")
 symbols(orstationc$lon, orstationc$lat, circles=orstationc$elev, inches=0.1, add=T)
 
-ggplot(orstationc, aes(x=lon, y=lat, size=elev)) + geom_point(shape=21, color="black", fill="lightblue")
+ggplot() + 
+  geom_sf(data = orcounty_sf, fill = NA) +
+  geom_point(data = orstationc, aes(x=lon, y=lat, size=elev), shape=21, color="black", fill="lightblue") +
+  coord_sf(xlim = c(-125, -116), ylim = c(41, 47), expand = FALSE) +
+  theme_bw()
 
 library(lattice)
 

R/plots_02.R

@@ -7,12 +7,87 @@ library(ncdf4)
 library(CFtime)
 library(RColorBrewer)
 
+library(lattice)
+
+attach(sumcr)
+coplot(WidthWS ~ DepthWS | Reach, pch=14+as.integer(Reach), cex=1.5, 
+       number=3, columns=3, overlap=0,
+       panel=function(x,y,...) {
+         panel.smooth(x,y,span=.8,iter=5,...)
+         abline(lm(y ~ x), col="blue")
+       } )
+
+ggplot(sumcr, aes(x=DepthWS, y=WidthWS)) +
+  stat_smooth(method = "lm") + 
+  geom_point() + 
+  facet_wrap(~ Reach) + 
+  theme(aspect.ratio = 1)
+
+detach(sumcr)
+
+# create some conditioning variables
+attach(yellpratio)
+Elevation <- equal.count(Elev,4,.25)
+Latitude <- equal.count(Lat,2,.25)
+Longitude <- equal.count(Lon,2,.25)
+
+# January vs July Precipitation Ratios by Elevation
+plot2 <- xyplot(APJan ~ APJul | Elevation,
+                layout = c(2, 2),
+                panel = function(x, y) {
+                  panel.grid(v=2)
+                  panel.xyplot(x, y)
+                  panel.abline(lm(y~x))
+                },
+                xlab = "APJul",
+                ylab = "APJan")
+print(plot2)
+
+# January vs July Precipitation Ratios by Latitude and Longitude
+plot3 <- xyplot(APJan ~ APJul | Latitude*Longitude, 
+                layout = c(2, 2),
+                panel = function(x, y) {
+                  panel.grid(v=2)
+                  panel.xyplot(x, y)
+                  panel.loess(x, y, span = .8, degree = 1, family="gaussian")
+                  panel.abline(lm(y~x))
+                },
+                xlab = "APJul",
+                ylab = "APJan")
+print(plot3)
+
+# create an elevation zones factor
+yellpratio$Elev_zones <- cut(Elevation, 4, labels=c("Elev1 (lowest)", "Elev2", "Elev3", "Elev4 (highest)"))
+
+ggplot(yellpratio, aes(x=APJul, y=APJan)) +
+  stat_smooth(method = "lm") + 
+  geom_point() + 
+  facet_wrap( ~ Elev_zones)
+
+# create conditioning variables
+Loc_factor <- rep(0, length(yellpratio$Lat))
+Loc_factor[(yellpratio$Lat >= median(yellpratio$Lat) & yellpratio$Lon < median(yellpratio$Lon))] <- 1
+Loc_factor[(yellpratio$Lat >= median(yellpratio$Lat) & yellpratio$Lon >= median(yellpratio$Lon))] <- 2
+Loc_factor[(yellpratio$Lat < median(yellpratio$Lat) & yellpratio$Lon < median(yellpratio$Lon))] <- 3
+Loc_factor[(yellpratio$Lat < median(yellpratio$Lat) & yellpratio$Lon >= median(yellpratio$Lon))] <- 4
+# convert to factor, and add level lables
+yellpratio$Loc_factor <- as.factor(Loc_factor)
+levels(yellpratio$Loc_factor) = c("Low Lon/High Lat", "High Lon/High Lat", "Low Lon/Low Lat", "High Lon/Low Lat")
+
+ggplot(yellpratio, aes(x=APJul, y=APJan)) +
+  stat_smooth(method = "loess", span = 0.9, col="red") + 
+  stat_smooth(method = "lm") + 
+  geom_point() + 
+  facet_wrap(~Loc_factor)
+
+
 # world_sf
-shapefile <- "/Users/bartlein/Dropbox/DataVis/working/data/shp_files/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp"
-world_sf <- st_read(shapefile)
+world_sf <- st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))
+world_sf
+plot(world_sf)
+
 world_otl_sf <- st_geometry(world_sf)
-plot(world_otl_sf) 
-class(world_otl_sf)
+plot(world_otl_sf)
 
 # ggplot map of world_outline
 ggplot() + 
@@ -37,8 +112,6 @@ ggplot()  +
   geom_stars(data = tree) +
   scale_fill_gradient(low = "white", high = "darkgreen", na.value = "transparent" ) +
   geom_sf(data = world_otl_sf, col = "black", fill = NA) +
-  # scale_x_continuous(breaks = breaks_x) +
-  # scale_y_continuous(breaks = breaks_y) +
   scale_x_continuous(breaks = seq(-180, 180, by = 30)) +
   scale_y_continuous(breaks = seq(-90, 90, by = 30)) +
   coord_sf(xlim = c(-180, +180), ylim = c(-90, 90), expand = FALSE) +
@@ -53,7 +126,6 @@ dname <- "tmp"  # note: tmp means temperature (not temporary)
 
 # stars read of netCDF
 tmp_stars <- read_ncdf(ncfname, var=dname)
-nc_close(ncin)
 
 # replace times
 #attr(tmp_stars, "dimensions")$time$values <- c(1:12)
@@ -62,14 +134,16 @@ attr(tmp_stars, "dimensions")$time$values <- c("Jan","Feb","Mar","Apr","May","Ju
 tmp_stars
 plot(tmp_stars)
 
+
+
 # ggplot2 map of tmp
 ggplot()  +
   geom_stars(data = tmp_stars) +
   scale_fill_gradient2(low = "darkblue", mid="white", high = "darkred", midpoint = 0) +
   geom_sf(data = world_otl_sf, col = "black", fill = NA) +
   facet_wrap("time", nrow = 4, ncol = 3) +
-  scale_x_continuous(breaks = breaks_x) +
-  scale_y_continuous(breaks = breaks_y) +
+  scale_x_continuous(breaks = seq(-180, 180, by = 30)) +
+  scale_y_continuous(breaks = seq(-90, 90, by = 30)) +
   coord_sf(xlim = c(-180, +175), ylim = c(-90, 90), expand = FALSE) +
   theme_bw() 
 
@@ -83,10 +157,12 @@ dim(tmp_array)
 nc_close(ncin)
 
 tmp_hovlat <- apply(tmp_array, c(2,3), mean, na.rm=TRUE)
-dim(tmp_hovlat); graphics::image(t(tmp_hovlat))
+dim(tmp_hovlat)
+image(t(tmp_hovlat))
 
 tmp_hovlon <- apply(tmp_array, c(1,3), mean, na.rm=TRUE)
-dim(tmp_hovlon); graphics::image(tmp_hovlon)
+dim(tmp_hovlon)
+image(tmp_hovlon)
 
 # Hovmoller plots
 
@@ -114,28 +190,12 @@ nlon * nlat
 cf <- CFtime(tunits$value, calendar = "proleptic_gregorian", time) # convert time to CFtime class
 cf
 timestamps <- CFtimestamp(cf) # get character-string times
-timestamps
-class(timestamps)
 time_cf <- CFparse(cf, timestamps) # parse the string into date components
 head(time_cf); tail(time_cf)
 class(time_cf)
 
-# check longitudes
-lon
-
-# rotate lons
-lontemp <- lon
-lon[1:(nlon/2)] <- lontemp[((nlon/2)+1):nlon] 
-lon[((nlon/2)+1):nlon] <- lontemp[1:(nlon/2)] + 360.0
-lon
-
-# rotate data while readling, get the west half of the array
-tmp_anm_array <- array(NA, c(nlon, nlat, nt))
-tmp_anm_array[1:(nlon/2),,] <- ncvar_get(ncin, dname, start = c(((nlon/2)+1), 1, 1), count = c(nlon/2, nlat, nt))
-# get the east half
-tmp_anm_array[((nlon/2)+1):nlon,,] <- ncvar_get(ncin, dname, start = c(1, 1, 1), count = c(nlon/2, nlat, nt))
-dim(tmp_anm_array)
-nc_close(ncin)
+# data
+tmp_anm_array <- ncvar_get(ncin, dname)
 
 # get begining year and ending year
 beg_yr <- time_cf$year[1]
@@ -167,12 +227,70 @@ ggplot()  +
   geom_tile(data = tmp_df01, aes(x = lon, y = lat, fill = tmp_anm)) +
   scale_fill_gradient2(low = "darkblue", mid="white", high = "darkred", midpoint = 0) +
   geom_sf(data = world_otl_sf, col = "black", fill = NA) +
-  scale_x_continuous(breaks = breaks_x) +
-  scale_y_continuous(breaks = breaks_y) +
-  coord_sf(xlim = c(-0, 360), ylim = c(-90, 90), expand = FALSE) +
+  scale_x_continuous(breaks = seq(-180, 180, by = 30)) +
+  scale_y_continuous(breaks = seq(-90, 90, by = 30)) +
+  coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
   labs(title=title, subtitle=pt1, fill="Anomalies") + 
   theme_bw() + theme(legend.position="bottom")
 
+####
+
+# rotate lons
+lon
+lontemp <- lon
+lon[1:(nlon/2)] <- lontemp[((nlon/2)+1):nlon] # 0 to + 180
+lon[((nlon/2)+1):nlon] <- lontemp[1:(nlon/2)] + 360.0 # 180 to 360
+lon
+
+# rotate data while reading...
+# get the west half of the new array
+tmp_anm_array <- array(NA, c(nlon, nlat, nt))
+tmp_anm_array[1:(nlon/2),,] <- ncvar_get(ncin, dname, start = c(((nlon/2)+1), 1, 1), count = c(nlon/2, nlat, nt))
+# get the east half
+tmp_anm_array[((nlon/2)+1):nlon,,] <- ncvar_get(ncin, dname, start = c(1, 1, 1), count = c(nlon/2, nlat, nt))
+dim(tmp_anm_array)
+
+# close the netCDF file
+nc_close(ncin)
+
+# vector of lons and lats
+lonlat <- as.matrix(expand.grid(lon,lat))
+dim(lonlat)
+head(lonlat); tail(lonlat)
+
+# vector of `tmp` values
+n <- 1957 # should be jan 2013
+tmp_vec <- as.vector(tmp_anm_array[,,n])
+length(tmp_vec)
+
+# create dataframe and add names
+tmp_df02 <- data.frame(cbind(lonlat,tmp_vec))
+names(tmp_df02) <- c("lon", "lat", "tmp_anm")
+
+library(maps)
+
+world2_sf <- st_as_sf(maps::map("world2", plot = FALSE, fill = TRUE))
+world2_sf
+plot(world2_sf)
+
+world2_otl_sf <- st_geometry(world2_sf)
+plot(world2_otl_sf)
+
+pt1 <- paste(as.character(time_cf$year[n]), as.character(time_cf$month[n]), sep = "-")
+title <- paste ("HadCRUTv5 -- 2m Air Temperature Anomalies")
+
+# ggplot2 map of tmp
+ggplot()  +
+  geom_tile(data = tmp_df02, aes(x = lon, y = lat, fill = tmp_anm)) +
+  scale_fill_gradient2(low = "darkblue", mid="white", high = "darkred", midpoint = 0) +
+  geom_sf(data = world2_otl_sf, col = "black", fill = NA) +
+  scale_x_continuous(breaks = seq(0, 360, by = 30)) +
+  scale_y_continuous(breaks = seq(-90, 90, by = 30)) +
+  coord_sf(xlim = c(0, 360), ylim = c(-90, 90), expand = FALSE) +
+  labs(title=title, subtitle=pt1, fill="Anomalies") + 
+  theme_bw() + theme(legend.position="bottom")
+
+
 # Hovmöller plots 
 tmp_hovlat <- apply(tmp_anm_array, c(2,3), mean, na.rm=TRUE)
 dim(tmp_hovlat)

docs/R/namer.R

@@ -1,5 +1,6 @@
 library(namer)
 
-filename <- "/Users/bartlein/Dropbox/DataVis/REarthSysSci/hdf5_intro.Rmd" # path to and name of .Rmd file 
+filename <- "/Users/bartlein/Dropbox/DataVis/REarthSysSci/plots1.Rmd" # path to and name of .Rmd file 
+unname_chunks(filename)
 name_chunks(filename)
 

docs/R/plots_01.R

@@ -9,7 +9,6 @@ library(ggplot2)
 attach(sumcr)
 plot(Length)
 
-qplot(seq(1:length(Length)), Length)
 ggplot(data=sumcr, aes(x = Length)) + 
   geom_point(aes(x = 1:nrow(sumcr), y = Length)) + 
   geom_line(aes(x = 1:nrow(sumcr), y = Length)) +
@@ -69,6 +68,7 @@ ggplot(specmap, aes(x = O18)) +
   geom_line(stat = "density") +
   geom_rug(data = specmap, aes(x = O18)) 
 
+
 # scatter plots
 # use Oregon climate-station data 
 
@@ -147,7 +147,11 @@ detach(sumcr)
 plot(orstationc$lon, orstationc$lat, type="n")
 symbols(orstationc$lon, orstationc$lat, circles=orstationc$elev, inches=0.1, add=T)
 
-ggplot(orstationc, aes(x=lon, y=lat, size=elev)) + geom_point(shape=21, color="black", fill="lightblue")
+ggplot() + 
+  geom_sf(data = orcounty_sf, fill = NA) +
+  geom_point(data = orstationc, aes(x=lon, y=lat, size=elev), shape=21, color="black", fill="lightblue") +
+  coord_sf(xlim = c(-125, -116), ylim = c(41, 47), expand = FALSE) +
+  theme_bw()
 
 library(lattice)