######################################################################################################## #------------------------------------------------------------------------------------------------------- # Author : A. John Woodill # Date : 1/26/2015 # Filename : California.R # Code : Precipitation Spatial Post for Website # Specifics : Merges mean precipitation for California plotted with Spatial object in ggplot2 #------------------------------------------------------------------------------------------------------- ######################################################################################################## library(plyr) library(dplyr) library(ggplot2) library(automap) # For Spacial Data library(akima) # For interpolation #------------------------------------------------------ # Data wrangling of PRISM data and Spatial Objects #------------------------------------------------------ CA_2000_appt <- read.csv("CA_2000_appt.csv") sub_data = as.data.frame(CA_2000_appt) coord_vars = c("latitude","longitude") data_vars = setdiff(colnames(sub_data), coord_vars) # All colnames except for Lat/Long sp_points = SpatialPoints(sub_data[,coord_vars]) # All lat/long in Spatial Point format sp_df = SpatialPointsDataFrame(sp_points, # Spatial Data Frame of lat/long and remaining cols sub_data[,data_vars,drop=FALSE]) regions <- c("california") map_base_data <- filter(map_data("state"), region == "california") # Get California lat/long coord. map_base_data <- rename(map_base_data, longitude = long) # Rename columns map_base_data <- rename(map_base_data, latitude = lat) # Createsa Spatial Polygon for each state with lat/long # ---> Used for plotting the map state_regions = function(x) { # Gets lat/long for each state to build Spatial Polygon state = unique(x$region) print(state) Polygons(list(Polygon(x[,c("latitude","longitude")])),ID=state) } state_pg = SpatialPolygons(dlply(map_base_data, # Builds a Spatial Polygon of all state/regions .(region), state_regions)) #------------------------------------------------------ # Spline interpolation with akima package #------------------------------------------------------ fld = with(sub_data, interp(x = longitude, y = latitude, z = APPT, duplicate="median", xo=seq(min(map_base_data$longitude), max(map_base_data$longitude), length = 100), yo=seq(min(map_base_data$latitude), max(map_base_data$latitude), length = 100), extrap=TRUE, linear=FALSE)) melt_x = rep(fld$x, times=length(fld$y)) # Pull out longitude values from list melt_y = rep(fld$y, each=length(fld$x)) # Pull out latitude values from list melt_z = as.vector(fld$z) # Pull out appt values from list level_data = data.frame(longitude=melt_x, # Build data.frame latitude=melt_y, APPT=melt_z) interp_data = na.omit(level_data) # Remove all NA values grid_points = SpatialPoints(interp_data[,2:1]) # Build Spatial Points into object in_points = !is.na(over(grid_points,state_pg)) # Logical determining points inside all regions and not NA inside_points = interp_data[in_points, ] # Removes all points outside of Spatial Polygons #------------------------------------------------------ # Plot Spatial Objects with ggplot2 #------------------------------------------------------ map_base_aesthetics = aes(x=longitude, y=latitude, group=group) # Aesthetics for ggplot2 map_base = geom_polygon(data=map_base_data, map_base_aesthetics) # Map Base borders = geom_polygon(data=map_base_data, # Draws boundaries map_base_aesthetics, color="black", fill=NA) ggplot(data=inside_points, aes(x=longitude, y=latitude)) + # Setup ggplot2 geom_tile(aes(fill=APPT)) + # Initial overlay of appt stat_contour(aes(z=APPT)) + # Create contours for concentrations coord_equal() + # Equalize plots scale_fill_gradient2(low="blue", mid="white",high="darkgreen", # Set colors for low, mid, high midpoint=mean(inside_points$APPT)) + ggtitle("California Precipitation - 2000") + # Plot title borders # Draw California border