El Paso Rain
##Getting Data from the Global Historical Climatology Website: https://www.ncdc.noaa.gov/data-access/land-based-station-data/land-based-datasets/global-historical-climatology-network-ghcn
#Load needed Libraries
suppressMessages(library(rnoaa))
suppressMessages(library(plyr))
suppressMessages(library(dplyr))
suppressMessages(library(raster))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))#State Boundaries To be used for mapping.
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)
TX = subset(States, ID == 2 | ID == 42 | ID == 35 | ID == 30)#World Boundaires Just for mapping.
World = map("world",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(World$names, ":"),
function(x) x[1])
World = map2SpatialPolygons(World, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(World, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(World),
row.names = pid)
World = SpatialPolygonsDataFrame(World, p.df)
World.Crop = crop(World, extent(TX))#Find Weather Stations Looking for weather stations within 250 km of El Paso. 56 Weather station are identified.
Stations = isd_stations_search(lat = 31.7619, lon = -106.4850, #El Paso Coordinates
radius = 250, bbox = NULL) #250km radius of El Paso Texas
Station.pnts = as.data.frame(Stations)
Station.pnts = SpatialPointsDataFrame(Station.pnts[,c("lon", "lat")], Station.pnts)
proj4string(Station.pnts) = proj4string(States)
#56 Weather station are identified.
Station.pnts## class : SpatialPointsDataFrame
## features : 56
## extent : -108.983, -104.258, 30.367, 33.9 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 12
## names : usaf, wban, station_name, ctry, state, icao, lat, lon, elev_m, begin, end, distance
## min values : 690014, 00345, ABRAHAM GONZALEZ INTL, MX, , , 30.367, -108.983, 985.1, 19410401, 19480701, 11.2709323917735
## max values : A00014, 99999, WHITE SANDS NORTH, US, TX, MMCS, 33.9, -104.258, 2810, 20120215, 20190925, 237.878450840898#Weather Station Location Map
Wrld_df = fortify(World.Crop)## Regions defined for each PolygonsTX_df = fortify(TX)## Regions defined for each PolygonsBr.plt = ggplot(Wrld_df,
aes(long,lat, group=group)) +
geom_polygon(fill="tan", col="black") +
geom_polygon(data = TX_df,
aes(long,lat, group=group),
fill="tan", col="black") +
geom_point(data=Station.pnts@data,
aes(lon, lat,
group=NULL),
col = "red", size=1, shape=2) +
xlab("Longitude") +
ylab("Latitude") +
coord_equal() +
ggtitle("Weather Stations within 250km") +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=22, face="bold", hjust = 0.5))
Br.plt
#Weather Data Pull Data for selected weather stations
El.Paso = data.frame(id = "El Paso", latitude = 31.7619, longitude = -106.4850) #Loaction of El Paso
#This code selcts the weather stations, date range, and variables wanted
station_data = ghcnd_stations()
Stat.df = meteo_nearby_stations(lat_lon_df = El.Paso,
station_data = station_data,
radius = 250, limit=250, #look for locations within 250km
var = "PRCP", #Precipitation data only (tenths of millimeters)
year_min = 2018, year_max = 2018) #Get data for year 2018
Station.info = as.data.frame(Stat.df) #record station identifying info
colnames(Station.info)[1] = "id"
Clim.df = meteo_pull_monitors(Stat.df$"El Paso"$id,
date_min = "2018-01-01",
date_max = "2018-12-31")
Clim.df = as.data.frame(Clim.df)
Keep = as.data.frame(Clim.df) %>%
filter(is.na(prcp)==FALSE) #Remove missing values
Keep = Keep[,c("id","date","prcp")]
Keep = plyr::join(Keep, Station.info, by = "id", type = "left")#Check and Save data Data file includes daily precipitation (prcp, in tenths of millimeters), staion id (id), Long, Lat and distance from El Paso.
head(Keep)## id date prcp El.Paso.name El.Paso.latitude
## 1 US1TXEP0053 1/11/2018 28 UTEP EHS 31.7661
## 2 US1TXEP0053 2/16/2018 25 UTEP EHS 31.7661
## 3 US1TXEP0053 2/17/2018 18 UTEP EHS 31.7661
## 4 US1TXEP0053 5/22/2018 61 UTEP EHS 31.7661
## 5 US1TXEP0053 6/4/2018 33 UTEP EHS 31.7661
## 6 US1TXEP0053 6/16/2018 99 UTEP EHS 31.7661
## El.Paso.longitude El.Paso.distance
## 1 -106.5048 1.929282
## 2 -106.5048 1.929282
## 3 -106.5048 1.929282
## 4 -106.5048 1.929282
## 5 -106.5048 1.929282
## 6 -106.5048 1.929282#Save data as a comma seperated text file
write.csv(Keep, "./Precip_093019.csv") #Save copy