1 |
f1c5b094
|
Benoit Parmentier
|
##################################### METHODS COMPARISON part 6 ##########################################
|
2 |
|
|
#################################### Spatial Analysis ############################################
|
3 |
|
|
#This script utilizes the R ojbects created during the interpolation phase. #
|
4 |
|
|
#At this stage the script produces figures of various accuracy metrics and compare methods: #
|
5 |
|
|
#This scripts focuses on a detailed studay of differences in the predictions of CAI_kr and FUsion_Kr #
|
6 |
|
|
#AUTHOR: Benoit Parmentier #
|
7 |
|
|
#DATE: 11/30/2012 #
|
8 |
|
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- #
|
9 |
|
|
###################################################################################################
|
10 |
|
|
|
11 |
|
|
###Loading R library and packages
|
12 |
|
|
library(gtools) # loading some useful tools such as mixedsort
|
13 |
|
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
14 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
15 |
|
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
16 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
17 |
|
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
18 |
|
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
19 |
|
|
library(spgwr)
|
20 |
|
|
library(gpclib)
|
21 |
|
|
library(maptools)
|
22 |
|
|
library(graphics)
|
23 |
|
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
24 |
|
|
library(raster)
|
25 |
|
|
library(rasterVis)
|
26 |
|
|
library(plotrix) #Draw circle on graph
|
27 |
|
|
library(reshape)
|
28 |
|
|
library(RCurl)
|
29 |
|
|
######### Functions used in the script
|
30 |
|
|
#loading R objects that might have similar names
|
31 |
|
|
|
32 |
|
|
download_files_with_pattern_ftp<-function(url_dir,out_path,file_pattern){
|
33 |
|
|
#sadfd
|
34 |
|
|
#adsfd
|
35 |
|
|
#library(RCurl) #modify later using require
|
36 |
|
|
# FTP
|
37 |
|
|
url<-url_dir
|
38 |
|
|
filenames = getURL(url, ftp.use.epsv = FALSE, dirlistonly = TRUE)
|
39 |
|
|
file_pattern_rx<-glob2rx(file_pattern)
|
40 |
|
|
# Deal with newlines as \n or \r\n. (BDR)
|
41 |
|
|
# Or alternatively, instruct libcurl to change \n's to \r\n's for us with crlf = TRUE
|
42 |
|
|
# filenames = getURL(url, ftp.use.epsv = FALSE, ftplistonly = TRUE, crlf = TRUE)
|
43 |
|
|
|
44 |
|
|
filenames = paste(url, strsplit(filenames, "\r*\n")[[1]], sep = "")
|
45 |
|
|
filenames_all<-filenames
|
46 |
|
|
#Now subset filenames to match wanted files...
|
47 |
|
|
i_file<-grep(file_pattern_rx,filenames_all,value=TRUE) # using grep with "value" extracts the matching names
|
48 |
|
|
pos<-match(i_file,filenames_all)
|
49 |
|
|
|
50 |
|
|
filenames<-filenames_all[pos]
|
51 |
|
|
53)}))
|
52 |
|
|
for(i in 1:length(filenames)){
|
53 |
|
|
out_file<-basename(filenames[i])
|
54 |
|
|
out_file_path<-file.path(out_path,out_file)
|
55 |
|
|
download.file(filenames[i],destfile=out_file_path)
|
56 |
|
|
#creating the file names
|
57 |
|
|
}
|
58 |
|
|
}
|
59 |
|
|
|
60 |
|
|
###Parameters and arguments
|
61 |
|
|
|
62 |
|
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
63 |
|
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
64 |
|
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
65 |
|
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
66 |
|
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
67 |
|
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
68 |
|
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
69 |
|
|
infile6<-"OR83M_state_outline.shp"
|
70 |
|
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
71 |
|
|
|
72 |
|
|
out_prefix<-"methods_11012012_"
|
73 |
|
|
|
74 |
|
|
path_wd<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012"
|
75 |
|
|
|
76 |
|
|
##### MAIN BODY OF CODE ######
|
77 |
|
|
|
78 |
|
|
setwd(path_wd)
|
79 |
|
|
latlong <- "+init=epsg:4326"
|
80 |
|
|
url_dir1<- "ftp://anonymous:test@ftp.wcc.nrcs.usda.gov/data/snow/snotel/snothist/"
|
81 |
|
|
url_dir2<- "ftp://anonymous:test@ftp.wcc.nrcs.usda.gov/data/snow/snotel/cards/oregon/"
|
82 |
|
|
file_pat<-"*or*.txt"
|
83 |
|
|
file_pat1<-"*listor.txt"
|
84 |
|
|
file_pat2<-"*all.txt"
|
85 |
|
|
output_path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012/snotel_data_11302012"
|
86 |
|
|
download_files_with_pattern_ftp(url_dir1,path_wd,file_pat1)
|
87 |
|
|
download_files_with_pattern_ftp(url_dir1,output_path,file_pat1)
|
88 |
|
|
download_files_with_pattern_ftp(url_dir2,output_path,file_pat2)
|
89 |
|
|
|
90 |
|
|
|
91 |
|
|
reg_outline<-readOGR(".",sub(".shp","",infile6))
|
92 |
|
|
loc_stat_file<-"listor.txt"
|
93 |
|
|
tmp<-read.table(loc_stat_file,sep=" ",skip=4)
|
94 |
|
|
|
95 |
|
|
colwidths<-c(5,3,4,5,9,9,5,6,6,50)
|
96 |
|
|
colab<- c("no", "ST", "CTY", "tyype", "HUC", "station", "Lat_pc", "Long_pc", "elev", "sitename")
|
97 |
|
|
loc_tab_snot<-read.fwf(file=loc_stat_file,widths=colwidths,skip=3)
|
98 |
|
|
names(loc_tab_snot)<-colab
|
99 |
|
|
|
100 |
|
|
y<-loc_stations_snot$Lat_pc/100
|
101 |
|
|
x<-(loc_stations_snot$Long_pc/100)*-1
|
102 |
|
|
|
103 |
|
|
tmp<-SpatialPointsDataFrame(coords=cbind(x,y),data=loc_stations_snot)
|
104 |
|
|
proj4string(tmp)<-CRS(latlong)
|
105 |
|
|
loc_stat_snot<-spTransform(tmp,CRS(proj4string(reg_outline)))
|
106 |
|
|
|
107 |
|
|
outfile_name<-"loc_stat_snot.shp"
|
108 |
|
|
writeOGR(obj=loc_stat_snot,layer="loc_stat_snot",dsn=outfile_name,driver="ESRI Shapefile", overwrite=T)
|
109 |
|
|
|
110 |
|
|
file_pat2_rx<-glob2rx(file_pat2)
|
111 |
|
|
snot_data_filename<-list.files(path=output_path,pattern=file_pat2_rx)
|
112 |
|
|
stat_id<-sub("_all.txt","",snot_data_filename)
|
113 |
|
|
tmpchr<-trim(levels(loc_stat_snot$station))
|
114 |
|
|
loc_stat_snot$station<-tolower(tmpchr)
|
115 |
|
|
|
116 |
|
|
m<-match(stat_id,loc_stat_snot$station) #This matches the station id location wtih the data files...
|
117 |
|
|
|
118 |
|
|
|
119 |
|
|
#put in a loop
|
120 |
|
|
#create giant table/data.frame from all the files...?
|
121 |
|
|
snot_data_name<-file.path(output_path,snot_data_filename[1])
|
122 |
|
|
snot_data_tab<-read.table(snot_data_name,sep="\t",skip=1) #This is a tab delimited file...
|
123 |
|
|
colab<- c("date", "pill", "prec","tmax", "tmin", "tavg", "prcp")
|
124 |
|
|
names(snot_data_tab)<-colab
|
125 |
|
|
|
126 |
|
|
snot_fnames<-file.path(output_path,snot_data_filename)
|
127 |
|
|
|
128 |
|
|
header=TRUE, sep="\t")}))
|
129 |
|
|
read_snot_data<-function(files){
|
130 |
|
|
snot_data_tab<-read.table(files,sep="\t",skip=1) #This is a tab delimited file...
|
131 |
|
|
stat_id_tmp<-sub("_all.txt","",files)
|
132 |
|
|
names(snot_data_tab)<-c("date", "pill", "prec", "tmax", "tmin", "tavg", "prcp")
|
133 |
|
|
snot_data_tab$id<-stat_id_tmp
|
134 |
|
|
return(snot_data_tab)
|
135 |
|
|
}
|
136 |
|
|
dataset<-do.call("rbind",lapply(file_list, FUN=read_snot_data))
|
137 |
|
|
|
138 |
|
|
dataset$id<-basename(dataset$id)
|
139 |
|
|
|
140 |
|
|
snot_loc_data_OR<-merge(loc_stat_snot,dataset,by.x="station",by.y="id")
|
141 |
|
|
|
142 |
|
|
tmp44<-subset(snot_loc_data_OR,station=="17d02s")
|
143 |
|
|
tail(snot_data_tab)
|
144 |
|
|
tmp45<-tmp44[tmp44$date=="93012",]
|
145 |
|
|
tail(tmp45)
|
146 |
|
|
tmp45<-tmp44[tmp44$date=="92912",]
|
147 |
|
|
|
148 |
|
|
test_data<-subset(snot_loc_data_OR,date=="90112" & station=="17d07s")
|
149 |
|
|
|
150 |
|
|
#Now subset for specific dates...
|
151 |
|
|
|
152 |
|
|
date_pat="*10" #anything in 2010...
|
153 |
|
|
date_pat_rx<-glob2rx(date_pat)
|
154 |
|
|
#grep(tmp_data$date
|
155 |
|
|
date_selected<-grep(date_pat_rx,snot_loc_data_OR$date,value=TRUE) # using grep with "value" extracts the matching names
|
156 |
|
|
#pos<-match(i_file,filenames_all)
|
157 |
|
|
|
158 |
|
|
snot_OR_2010<-subset(snot_loc_data_OR,date%in%date_selected)
|
159 |
|
|
|
160 |
|
|
#tansform in SPDF
|
161 |
|
|
coords<-cbind(snot_OR_2010$x,snot_OR_2010$y)
|
162 |
|
|
snot_OR_2010_sp<-SpatialPointsDataFrame(coords=coords,data=snot_OR_2010)
|
163 |
|
|
length(unique(snot_OR_2010_sp$station)) #Find out about the number of stations...
|
164 |
|
|
dim(snot_OR_2010_sp)
|
165 |
|
|
#write out in SPDF
|
166 |
|
|
outfile_name<-paste("snot_OR_2010_sp2_",out_prefix,".shp",sep="")
|
167 |
|
|
writeOGR(obj=snot_OR_2010_sp,layer="snot_OR_2010",dsn=outfile_name,driver="ESRI Shapefile", overwrite=T)
|
168 |
|
|
|
169 |
|
|
#####END OF SCRIPT
|
170 |
|
|
#
|
171 |
|
|
# dataset <- do.call("rbind",lapply(file_list,
|
172 |
|
|
# FUN=function(files){read.table(files,
|
173 |
|
|
# skip=1, sep="\t")}))
|
174 |
|
|
#
|
175 |
|
|
|
176 |
|
|
|