1 |
147da66d
|
Adam M. Wilson
|
###################################################################################
|
2 |
|
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
3 |
|
|
|
4 |
|
|
|
5 |
|
|
## connect to server of choice
|
6 |
|
|
#system("ssh litoria")
|
7 |
|
|
#R
|
8 |
|
|
|
9 |
|
|
library(sp)
|
10 |
|
|
library(rgdal)
|
11 |
|
|
library(reshape)
|
12 |
|
|
library(ncdf4)
|
13 |
|
|
library(geosphere)
|
14 |
|
|
library(rgeos)
|
15 |
|
|
library(multicore)
|
16 |
|
|
library(raster)
|
17 |
|
|
library(lattice)
|
18 |
|
|
library(rgl)
|
19 |
|
|
library(hdf5)
|
20 |
|
|
|
21 |
|
|
X11.options(type="Xlib")
|
22 |
|
|
ncores=20 #number of threads to use
|
23 |
|
|
|
24 |
|
|
setwd("/home/adamw/personal/projects/interp")
|
25 |
|
|
setwd("/home/adamw/acrobates/projects/interp")
|
26 |
|
|
|
27 |
|
|
roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")
|
28 |
|
|
roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
29 |
|
|
|
30 |
|
|
### Downloading data from LAADSWEB
|
31 |
|
|
# subset by geographic area of interest
|
32 |
|
|
# subset: 40-47, -115--125
|
33 |
|
|
|
34 |
|
|
## download data from ftp site. Unfortunately the selection has to be selected using the website and orders downloaded via ftp.
|
35 |
|
|
|
36 |
|
|
system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")
|
37 |
|
|
|
38 |
|
|
|
39 |
|
|
gdir="output/"
|
40 |
|
|
datadir="data/modis/MOD06_L2_hdf"
|
41 |
|
|
outdir="data/modis/MOD06_L2_nc"
|
42 |
|
|
|
43 |
|
|
fs=data.frame(
|
44 |
|
|
path=list.files(datadir,full=T,recursive=T,pattern="hdf"),
|
45 |
|
|
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))
|
46 |
|
|
fs$date=as.Date(substr(fs$file,11,17),"%Y%j")
|
47 |
|
|
fs$time=substr(fs$file,19,22)
|
48 |
|
|
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))
|
49 |
|
|
fs$path=as.character(fs$path)
|
50 |
|
|
fs$file=as.character(fs$file)
|
51 |
|
|
|
52 |
|
|
## output ROI
|
53 |
|
|
#get bounding box of region in m
|
54 |
|
|
ge=SpatialPoints(data.frame(lon=c(-125,-115),lat=c(40,47)))
|
55 |
|
|
projection(ge)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
56 |
|
|
ge2=spTransform(ge, CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
57 |
|
|
|
58 |
|
|
|
59 |
|
|
## vars
|
60 |
|
|
vars=paste(c(
|
61 |
|
|
"Cloud_Effective_Radius",
|
62 |
|
|
"Cloud_Effective_Radius_Uncertainty",
|
63 |
|
|
"Cloud_Optical_Thickness",
|
64 |
|
|
"Cloud_Optical_Thickness_Uncertainty",
|
65 |
|
|
"Cloud_Water_Path",
|
66 |
|
|
"Cloud_Water_Path_Uncertainty",
|
67 |
|
|
"Cloud_Phase_Optical_Properties",
|
68 |
|
|
"Cloud_Multi_Layer_Flag",
|
69 |
|
|
"Cloud_Mask_1km",
|
70 |
|
|
"Quality_Assurance_1km"))
|
71 |
|
|
|
72 |
|
|
#vars=paste(c("Cloud_Effective_Radius","Cloud_Optical_Thickness","Cloud_Water_Path"))
|
73 |
|
|
|
74 |
|
|
### Installation of hegtool
|
75 |
|
|
## needed 32-bit libraries and java for program to install correctly
|
76 |
|
|
|
77 |
|
|
# system(paste("hegtool -h ",fs$path[1],sep=""))
|
78 |
|
|
|
79 |
|
|
|
80 |
|
|
#### Function to generate hegtool parameter file for multi-band HDF-EOS file
|
81 |
|
|
swath2grid=function(i=1,files,vars=vars,outdir,upleft="47 -125",lowright="41 -115"){
|
82 |
|
|
file=fs$path[i]
|
83 |
|
|
print(paste("Starting file",basename(file)))
|
84 |
|
|
tempfile=paste(tempdir(),"/",sub("hdf$","tif",fs$file[i]),sep="")
|
85 |
|
|
date=fs$date[1]
|
86 |
|
|
origin=as.POSIXct("1970-01-01 00:00:00",tz="GMT")
|
87 |
|
|
### First write the parameter file (careful, heg is very finicky!)
|
88 |
|
|
hdr=paste("NUM_RUNS = ",length(vars),"|MULTI_BAND_GEOTIFF:",length(vars),sep="")
|
89 |
|
|
# hdr=paste("NUM_RUNS = ",length(vars),"|",sep="")
|
90 |
|
|
grp=paste("
|
91 |
|
|
BEGIN
|
92 |
|
|
INPUT_FILENAME=",file,"
|
93 |
|
|
OBJECT_NAME=mod06
|
94 |
|
|
FIELD_NAME=",vars,"|
|
95 |
|
|
BAND_NUMBER = 1
|
96 |
|
|
OUTPUT_PIXEL_SIZE_X=1000
|
97 |
|
|
OUTPUT_PIXEL_SIZE_Y=1000
|
98 |
|
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
99 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
100 |
|
|
RESAMPLING_TYPE = ",ifelse(grepl("Flag|Mask",vars),"NN","NN"),"
|
101 |
|
|
OUTPUT_PROJECTION_TYPE = SIN
|
102 |
|
|
ELLIPSOID_CODE = WGS84
|
103 |
|
|
OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
|
104 |
|
|
OUTPUT_TYPE = GEO
|
105 |
|
|
#OUTPUT_FILENAME = ",paste(tempfile,paste(vars,".tif",sep=""),sep=""),"
|
106 |
|
|
OUTPUT_FILENAME = ",tempfile,"
|
107 |
|
|
END
|
108 |
|
|
|
109 |
|
|
|
110 |
|
|
",sep="")
|
111 |
|
|
## if any remnants from previous runs remain, delete them
|
112 |
|
|
if(length(list.files(tempdir(),pattern=basename(tempfile)))>0)
|
113 |
|
|
file.remove(list.files(tempdir(),pattern=basename(tempfile),full=T))
|
114 |
|
|
## write it to a file
|
115 |
|
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
116 |
|
|
## now run the swath2grid tool - must be run as root (argh!)!
|
117 |
|
|
## write the tiff file
|
118 |
|
|
log=system(paste("sudo MRTDATADIR=/usr/local/heg/data ",
|
119 |
|
|
"PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",
|
120 |
|
|
paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
121 |
|
|
## convert to netcdf to extract metadata
|
122 |
|
|
system(paste("NCARG_ROOT=\"/usr/local/ncl_ncarg-6.0.0\" ncl_convert2nc ",
|
123 |
|
|
file," -o ",tempdir()," -B ",sep=""))
|
124 |
|
|
nc=nc_open(sub("[.]tif$",".nc",tempfile))
|
125 |
|
|
|
126 |
|
|
## read variables one by one into R to expand the file
|
127 |
|
|
for(v in vars){
|
128 |
|
|
atts=ncatt_get(nc,v)
|
129 |
|
|
outfile=sub("[.]hdf$",paste("_",v,".nc",sep=""),paste(outdir,"/",fs$file[i],sep=""))
|
130 |
|
|
QA=grepl("Flag|Mask",v)
|
131 |
|
|
td=expand(raster(tempfile,band=which(v==vars)),ge2,value=ifelse(QA,0,-9999))
|
132 |
|
|
projection(td)=CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0")
|
133 |
|
|
|
134 |
|
|
NAvalue(td)=ifelse(QA,0,-9999)
|
135 |
|
|
## get netcdf attributes
|
136 |
|
|
ncfile1=paste(tempdir(),"/",sub("hdf$","nc",basename(file)),sep="")
|
137 |
|
|
writeRaster(td,ncfile1,overwrite=T,datatype=ifelse(QA,"INT1U","INT2S"),NAflag=ifelse(QA,0,-9999))
|
138 |
|
|
|
139 |
|
|
## add time variable
|
140 |
|
|
print("Adding time dimension")
|
141 |
|
|
system(paste("ncecat -O -u time ",ncfile1,outfile))
|
142 |
|
|
system(paste("ncap2 -s \'time[time]=",
|
143 |
|
|
as.numeric(difftime(fs$datetime[i],origin,units="mins")),"\' ",outfile,sep=""))
|
144 |
|
|
|
145 |
|
|
|
146 |
|
|
######################################################################################
|
147 |
|
|
|
148 |
|
|
print("Updating netCDF dimensions")
|
149 |
|
|
## rename dimension variables
|
150 |
|
|
## writeraster incorrectly labels the dimensions! y=easting!
|
151 |
|
|
system(paste("ncrename -d northing,x -d easting,y -v northing,x -v easting,y -v variable,",v," ",outfile,sep=""))
|
152 |
|
|
system(paste("ncap2 -s \'lat[x,y]=0;lon[x,y]=0\' ",outfile))
|
153 |
|
|
system(paste("ncap2 -s \'sinusoidal=0\' ",outfile))
|
154 |
|
|
|
155 |
|
|
## Get corner coordinates and convert to cell centers
|
156 |
|
|
ncd=system(paste("gdalinfo ",tempfile),intern=T)
|
157 |
|
|
# UL= as.numeric(strsplit(gsub("[a-z]|[A-Z]|[\\]|[\t]|[\"]|[=]|[(]|[)]|[,]","",ncd[grep("Upper Left",ncd)]),
|
158 |
|
|
# " +")[[1]][2:3])+c(500,-500)
|
159 |
|
|
# LR= as.numeric(strsplit(gsub("[a-z]|[A-Z]|[\\]|[\t]|[\"]|[=]|[(]|[)]|[,]","",ncd[grep("Lower Right",ncd)]),
|
160 |
|
|
# " +")[[1]][2:3])+c(-500,500)
|
161 |
|
|
UL=c(extent(td)@xmin,extent(td)@ymax)+c(500,-500)
|
162 |
|
|
LR=c(extent(td)@xmax,extent(td)@ymin)+c(-500,500)
|
163 |
|
|
xvar=seq(UL[1],LR[1],by=1000)
|
164 |
|
|
yvar=seq(UL[2],LR[2],by=-1000)
|
165 |
|
|
nc=nc_open(outfile,write=T)
|
166 |
|
|
ncvar_put(nc,"x",vals=xvar)
|
167 |
|
|
ncvar_put(nc,"y",vals=yvar)
|
168 |
|
|
|
169 |
|
|
|
170 |
|
|
## add lat-lon grid
|
171 |
|
|
grid=expand.grid(x=xvar,y=yvar)
|
172 |
|
|
coordinates(grid)=c("x","y")
|
173 |
|
|
projection(grid)=CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0")
|
174 |
|
|
grid2=spTransform(grid,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
175 |
|
|
grid=SpatialPointsDataFrame(grid,data=cbind.data.frame(coordinates(grid),coordinates(grid2)))
|
176 |
|
|
gridded(grid)=T
|
177 |
|
|
fullgrid(grid)=T
|
178 |
|
|
colnames(grid@data)=c("x","y","lon","lat")
|
179 |
|
|
|
180 |
|
|
xlon=as.matrix(grid["lon"])
|
181 |
|
|
ylat=as.matrix(grid["lat"])
|
182 |
|
|
|
183 |
|
|
ncvar_put(nc,"lon",vals=xlon)
|
184 |
|
|
ncvar_put(nc,"lat",vals=ylat)
|
185 |
|
|
|
186 |
|
|
nc_close(nc)
|
187 |
|
|
|
188 |
|
|
print("Updating attributes")
|
189 |
|
|
system(paste("ncatted -a coordinates,",v,",o,c,\"lat lon\" ",outfile,sep=""))
|
190 |
|
|
system(paste("ncatted -a grid_mapping,",v,",o,c,\"sinusoidal\" ",outfile,sep=""))
|
191 |
|
|
|
192 |
|
|
system(paste("ncatted -a units,time,o,c,\"Minutes since ",origin,"\" ",outfile))
|
193 |
|
|
system(paste("ncatted -a long_name,time,o,c,\"time\"",outfile))
|
194 |
|
|
|
195 |
|
|
system(paste("ncatted -a units,y,o,c,\"m\" ",outfile))
|
196 |
|
|
system(paste("ncatted -a long_name,y,o,c,\"y coordinate of projection\" ",outfile))
|
197 |
|
|
system(paste("ncatted -a standard_name,y,o,c,\"projection_y_coordinate\" ",outfile))
|
198 |
|
|
|
199 |
|
|
system(paste("ncatted -a units,x,o,c,\"m\" ",outfile))
|
200 |
|
|
system(paste("ncatted -a long_name,x,o,c,\"x coordinate of projection\" ",outfile))
|
201 |
|
|
system(paste("ncatted -a standard_name,x,o,c,\"projection_x_coordinate\" ",outfile))
|
202 |
|
|
|
203 |
|
|
## grid attributes
|
204 |
|
|
system(paste("ncatted -a units,lat,o,c,\"degrees_north\" ",outfile))
|
205 |
|
|
system(paste("ncatted -a long_name,lat,o,c,\"latitude coordinate\" ",outfile))
|
206 |
|
|
system(paste("ncatted -a standard_name,lat,o,c,\"latitude\" ",outfile))
|
207 |
|
|
|
208 |
|
|
system(paste("ncatted -a units,lon,o,c,\"degrees_east\" ",outfile))
|
209 |
|
|
system(paste("ncatted -a long_name,lon,o,c,\"longitude coordinate\" ",outfile))
|
210 |
|
|
system(paste("ncatted -a standard_name,lon,o,c,\"longitude\" ",outfile))
|
211 |
|
|
|
212 |
|
|
system(paste("ncatted -a grid_mapping_name,sinusoidal,o,c,\"sinusoidal\" ",outfile))
|
213 |
|
|
system(paste("ncatted -a standard_parallel,sinusoidal,o,c,\"0\" ",outfile))
|
214 |
|
|
system(paste("ncatted -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",outfile))
|
215 |
|
|
system(paste("ncatted -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",outfile))
|
216 |
|
|
|
217 |
|
|
|
218 |
|
|
print(paste("Finished ", file))
|
219 |
|
|
|
220 |
|
|
}
|
221 |
|
|
## clean up and delete temporary files
|
222 |
|
|
if(length(list.files(tempdir(),pattern=basename(tempfile)))>0)
|
223 |
|
|
file.remove(list.files(tempdir(),pattern=basename(tempfile),full=T))
|
224 |
|
|
|
225 |
|
|
}
|
226 |
|
|
|
227 |
|
|
|
228 |
|
|
### update fs with completed files
|
229 |
|
|
fs$complete=fs$file%in%sub("nc$","hdf",list.files("data/modis/MOD06_L2_nc",pattern="nc$"))
|
230 |
|
|
table(fs$complete)
|
231 |
|
|
|
232 |
|
|
#### Run the gridding procedure
|
233 |
|
|
|
234 |
|
|
#for(fi in which(!fs$complete))
|
235 |
|
|
#swath2grid(fi,vars=vars,files=fs,outdir=outdir,upleft="47 -125",lowright="40 -115")
|
236 |
|
|
|
237 |
|
|
system("sudo ls")
|
238 |
|
|
|
239 |
|
|
mclapply(which(!fs$complete)[1:10],function(fi){
|
240 |
|
|
swath2grid(fi,vars=vars,files=fs,
|
241 |
|
|
outdir="data/modis/MOD06_L2_test",
|
242 |
|
|
upleft="47 -125",lowright="40 -115")},
|
243 |
|
|
mc.cores=20)
|
244 |
|
|
|
245 |
|
|
|
246 |
|
|
|
247 |
|
|
#### Merge the files
|
248 |
|
|
system(paste("cdo mergetime ",paste(list.files(outdir,full=T),collapse=" "),"data/modis/MOD06L2.nc"))
|
249 |
|
|
|
250 |
|
|
|
251 |
|
|
|
252 |
|
|
|
253 |
|
|
|
254 |
|
|
### subset to particular variables and drop uncertainties
|
255 |
|
|
f=fs[grep(vars[3],fs$file),]
|
256 |
|
|
f=f[!grepl("Uncertainty",f$file),]
|
257 |
|
|
|
258 |
|
|
getextent=function(file=f$file[1]) {
|
259 |
|
|
ext=extent(raster(file))
|
260 |
|
|
data.frame(file=file,xmin=attr(ext,"xmin"),xmax=attr(ext,"xmax"),ymin=attr(ext,"ymin"),ymax=attr(ext,"ymax"))
|
261 |
|
|
}
|
262 |
|
|
|
263 |
|
|
### some extents don't line up, find which ones...
|
264 |
|
|
ext=do.call(rbind.data.frame,mclapply(f$path,getextent))
|
265 |
|
|
f[which(ext$xmax!=ext$xmax[1]),]
|
266 |
|
|
f[which(ext$ymin!=ext$ymin[1]),]
|
267 |
|
|
|
268 |
|
|
### get extent of first image
|
269 |
|
|
ext1=extent(raster(f$path[1]))
|
270 |
|
|
### generate brick of all images
|
271 |
|
|
d=brick(lapply(f$path, function(i) {print(i) crop(raster(i),ext1)}))
|
272 |
|
|
|
273 |
|
|
|
274 |
|
|
## rescale data and identify NAs
|
275 |
|
|
NAvalue(d)=-9999
|
276 |
|
|
d=d*0.009999999776482582
|
277 |
|
|
|
278 |
|
|
plot(d)
|
279 |
|
|
|
280 |
|
|
## calculate overall mean
|
281 |
|
|
dm=mean(d,na.rm=T)
|
282 |
|
|
|
283 |
|
|
## plot it
|
284 |
|
|
X11.options(type="Xlib")
|
285 |
|
|
|
286 |
|
|
pdf("output/MOD06_mean.pdf",width=1000,height=600)
|
287 |
|
|
plot(d)
|
288 |
|
|
plot(roi,add=T)
|
289 |
|
|
dev.off()
|
290 |
|
|
|
291 |
|
|
|
292 |
|
|
|
293 |
|
|
|
294 |
|
|
|
295 |
|
|
|
296 |
|
|
|
297 |
|
|
|
298 |
|
|
|
299 |
|
|
|
300 |
|
|
|
301 |
|
|
|
302 |
|
|
|
303 |
|
|
#### load the functions
|
304 |
|
|
source("code/GHCN_functions.r")
|
305 |
|
|
|
306 |
|
|
|
307 |
|
|
|
308 |
|
|
|
309 |
|
|
|
310 |
|
|
#################################################
|
311 |
|
|
#### Download Data
|
312 |
|
|
dir.create("data/lst")
|
313 |
|
|
ftpsite="atlas.nceas.ucsb.edu:/home/parmentier/data_Oregon_stations/"
|
314 |
|
|
system(paste("scp -r wilson@",ftpsite,"mean_month* data/lst ",sep=""))
|
315 |
|
|
|
316 |
|
|
lst=brick(as.list(list.files("data/lst/",pattern=".*rescaled[.]rst",full=T)[c(4:12,1:3)]))
|
317 |
|
|
projection(lst)=CRS("+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")
|
318 |
|
|
lst=lst-272.15
|
319 |
|
|
|
320 |
|
|
|
321 |
|
|
#################################################
|
322 |
|
|
### Load PRISM
|
323 |
|
|
prism=brick("data/prism/prism_tmax.nc")
|
324 |
|
|
prism=projectRaster(prism,lst)
|
325 |
|
|
|
326 |
|
|
prism2=as(prism,"SpatialGridDataFrame");colnames(prism2@data)=paste("X",1:12,sep="")
|
327 |
|
|
prism2$source="prism"
|
328 |
|
|
prism2@data[c("lon","lat")]=coordinates(prism2)
|
329 |
|
|
|
330 |
|
|
lst2=as(lst,"SpatialGridDataFrame");colnames(lst2@data)=paste("X",1:12,sep="")
|
331 |
|
|
lst2$source="modis"
|
332 |
|
|
lst2@data[c("lon","lat")]=coordinates(lst2)
|
333 |
|
|
|
334 |
|
|
d=rbind.data.frame(melt(prism2@data,id.vars=c("source","lat","lon")),melt(lst2@data,id.vars=c("source","lat","lon")))
|
335 |
|
|
d$source=as.factor(d$source)
|
336 |
|
|
colnames(d)[grep("variable",colnames(d))]="month"
|
337 |
|
|
levels(d$month)=month.name
|
338 |
|
|
d=d[!is.na(d$value),]
|
339 |
|
|
|
340 |
|
|
d2=cast(d,lat+lon+month~source); gc()
|
341 |
|
|
d2$modis[d2$modis<(-50)]=NA
|
342 |
|
|
|
343 |
|
|
|
344 |
|
|
|
345 |
|
|
|
346 |
|
|
plot(subset(lst,1),col=rev(heat.colors(20)))
|
347 |
|
|
|
348 |
|
|
|
349 |
|
|
|
350 |
|
|
#### Explore prism-LST relationship (and land cover!?!)
|
351 |
|
|
|
352 |
|
|
png(width=1024,height=768,file="LSTvsPRISM.png")
|
353 |
|
|
xyplot(modis~prism|month,data=d2,panel=function(x,y,subscripts){
|
354 |
|
|
panel.xyplot(x,y,cex=.2,pch=16,col="black")
|
355 |
|
|
lm1=lm(y~x)
|
356 |
|
|
panel.abline(lm1)
|
357 |
|
|
panel.abline(0,1,col="red")
|
358 |
|
|
panel.text(-0,40,paste("R2=",round(summary(lm1)$r.squared,2)))
|
359 |
|
|
},ylab="MODIS Daytime LST (C)",xlab="PRISM Maximum Temperature (C)",
|
360 |
|
|
sub="Red line is y=x, black line is regression")
|
361 |
|
|
dev.off()
|
362 |
|
|
|
363 |
|
|
|
364 |
|
|
### Bunch of junk below here!
|
365 |
|
|
|
366 |
|
|
|
367 |
|
|
#### Perspective plot
|
368 |
|
|
open3d()
|
369 |
|
|
|
370 |
|
|
|
371 |
|
|
### load data to show some points
|
372 |
|
|
load("stroi.Rdata")
|
373 |
|
|
load("data/ghcn/roi_ghcn.Rdata")
|
374 |
|
|
d2=d[d$date=="2010-01-01"&d$variable=="tmax",]
|
375 |
|
|
|
376 |
|
|
r=extent(212836,234693,320614,367250)
|
377 |
|
|
z=as.matrix(raster(crop(lst,r),layer=1))
|
378 |
|
|
|
379 |
|
|
x <- 1:nrow(z)
|
380 |
|
|
y <- 1:ncol(z)
|
381 |
|
|
|
382 |
|
|
|
383 |
|
|
ncol=50
|
384 |
|
|
colorlut <- heat.colors(ncol,alpha=0) # height color lookup table
|
385 |
|
|
#col <- colorlut[ z-min(z)+1 ] # assign colors to heights for each point
|
386 |
|
|
col <- colorlut[trunc((z/max(z))*(ncol-1))+1] # assign colors to heights for each point
|
387 |
|
|
|
388 |
|
|
|
389 |
|
|
rgl.surface(x, y, z, color=col, alpha=0.75, back="lines")
|
390 |
|
|
rgl.points()
|
391 |
|
|
|
392 |
|
|
#### 2d example
|
393 |
|
|
n=100
|
394 |
|
|
xlim=c(-2,3)
|
395 |
|
|
x=seq(xlim[1],xlim[2],length=n)
|
396 |
|
|
e=rnorm(n,0,.5)
|
397 |
|
|
|
398 |
|
|
### climate function
|
399 |
|
|
yclimate<-function(x) x^4-x^3-7*x^2+1.2*x+15+(1.2*sin(x*.5))+(1*cos(x*5))
|
400 |
|
|
### daily weather function
|
401 |
|
|
yweather<-function(x) yclimate(x)+2.5*x+.3*x^3
|
402 |
|
|
|
403 |
|
|
y1=yclimate(x)
|
404 |
|
|
y2=yweather(x)
|
405 |
|
|
|
406 |
|
|
#points
|
407 |
|
|
s=c(-1.5,-.2,1.4,2.5)
|
408 |
|
|
s1=yclimate(s)
|
409 |
|
|
s2=yweather(s)
|
410 |
|
|
|
411 |
|
|
png(width=1024,height=768,file="anomalyapproach_%d.png",pointsize=28)
|
412 |
|
|
for(i in 1:6){
|
413 |
|
|
par(mfrow=c(2,1),mar=c(0,5,4,1))
|
414 |
|
|
plot(y2~x,type="l",xaxt="n",xlab="Space",ylab="Temperature",las=1,ylim=c(-9,22),
|
415 |
|
|
yaxp=c(0,20,1),col="transparent")
|
416 |
|
|
if(i<5) mtext("Space",1)
|
417 |
|
|
## points
|
418 |
|
|
if(i>=1) {
|
419 |
|
|
points(s,s2,pch=16,col="blue")
|
420 |
|
|
text(0.1,11,"Station data",col="blue")
|
421 |
|
|
}
|
422 |
|
|
if(i>=2) {
|
423 |
|
|
lines(y2~x,lwd=2)
|
424 |
|
|
points(s,s2,pch=16,col="blue") #put points back on top
|
425 |
|
|
text(xlim[1]+.5,-3,"Weather")
|
426 |
|
|
}
|
427 |
|
|
if(i>=3) {
|
428 |
|
|
lines(y1~x,col="darkgreen",lwd=2)
|
429 |
|
|
text(xlim[1]+.5,10,"Climate",col="darkgreen")
|
430 |
|
|
}
|
431 |
|
|
if(i>=4) segments(s,s2,s,s1,col="blue",lwd=2)
|
432 |
|
|
### anomalies
|
433 |
|
|
if(i>=5) {
|
434 |
|
|
par(mar=c(4,5,1,1))
|
435 |
|
|
plot(I(y2-y1)~x,type="l",xaxt="n",xlab="Space",
|
436 |
|
|
ylab="Temperature Anomaly \n (Weather-Climate)",las=1,ylim=c(-9,20),
|
437 |
|
|
yaxp=c(0,20,1),col="transparent")
|
438 |
|
|
## points
|
439 |
|
|
points(s,s2-s1,pch=16,col="blue")
|
440 |
|
|
abline(h=0,col="grey",lwd=2)
|
441 |
|
|
segments(s,0,s,s2-s1,col="blue",lwd=2)
|
442 |
|
|
}
|
443 |
|
|
if(i>=6) lines(I(y2-y1)~x,lwd=2)
|
444 |
|
|
}
|
445 |
|
|
dev.off()
|
446 |
|
|
|
447 |
|
|
|
448 |
|
|
|
449 |
|
|
##########################################################
|
450 |
|
|
#### Identify region of interest
|
451 |
|
|
|
452 |
|
|
### develop in region of interest spatial polygon
|
453 |
|
|
roi=readOGR("data/boundaries/statesp020.shp","statesp020")
|
454 |
|
|
proj4string(roi)=CRS("+proj=longlat +datum=WGS84")
|
455 |
|
|
roi=roi[roi$STATE=="Oregon",]
|
456 |
|
|
## buffer region of interest to include surrounding stations (in km)
|
457 |
|
|
roib=bufferroi(roi,distance=100)
|