Project

General

Profile

« Previous | Next » 

Revision f2b2e7d5

Added by Adam M. Wilson over 12 years ago

Added initial version of LST_Landcover exploration (not finished). Also moved prior version of interpolation procedure (bayesian krig) into repository.

View differences:

climate/extra/Interpolation_flowchart/flowchart.svg
1
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
2
<!-- Created with Inkscape (http://www.inkscape.org/) -->
3

  
4
<svg
5
   xmlns:dc="http://purl.org/dc/elements/1.1/"
6
   xmlns:cc="http://creativecommons.org/ns#"
7
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
8
   xmlns:svg="http://www.w3.org/2000/svg"
9
   xmlns="http://www.w3.org/2000/svg"
10
   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
11
   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
12
   width="633.65942"
13
   height="732.85431"
14
   id="svg2"
15
   version="1.1"
16
   inkscape:version="0.48.3.1 r9886"
17
   sodipodi:docname="flowchart.svg"
18
   inkscape:export-filename="/home/adamw/work/acrobates/data/personal/adamw/projects/environmental-layers/climate/extra/Interpolation_flowchart/flowchart.png"
19
   inkscape:export-xdpi="150"
20
   inkscape:export-ydpi="150">
21
  <sodipodi:namedview
22
     id="base"
23
     pagecolor="#ffffff"
24
     bordercolor="#666666"
25
     borderopacity="1.0"
26
     inkscape:pageopacity="1"
27
     inkscape:pageshadow="2"
28
     inkscape:zoom="0.98994949"
29
     inkscape:cx="349.46847"
30
     inkscape:cy="378.36894"
31
     inkscape:document-units="px"
32
     inkscape:current-layer="layer1"
33
     showgrid="false"
34
     inkscape:window-width="1280"
35
     inkscape:window-height="999"
36
     inkscape:window-x="0"
37
     inkscape:window-y="0"
38
     inkscape:window-maximized="1"
39
     fit-margin-top="0"
40
     fit-margin-left="0"
41
     fit-margin-right="0"
42
     fit-margin-bottom="0" />
43
  <defs
44
     id="defs4">
45
    <marker
46
       inkscape:stockid="Arrow1Lend"
47
       orient="auto"
48
       refY="0"
49
       refX="0"
50
       id="Arrow1Lend"
51
       style="overflow:visible">
52
      <path
53
         id="path4557"
54
         d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
55
         style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
56
         transform="matrix(-0.8,0,0,-0.8,-10,0)"
57
         inkscape:connector-curvature="0" />
58
    </marker>
59
    <marker
60
       inkscape:stockid="Arrow1LendF"
61
       orient="auto"
62
       refY="0"
63
       refX="0"
64
       id="Arrow1LendF"
65
       style="overflow:visible">
66
      <path
67
         id="path5295"
68
         d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
69
         style="fill:#818181;fill-rule:evenodd;stroke:#818181;stroke-width:1pt"
70
         transform="matrix(-0.8,0,0,-0.8,-10,0)"
71
         inkscape:connector-curvature="0" />
72
    </marker>
73
    <marker
74
       inkscape:stockid="Arrow1LendI"
75
       orient="auto"
76
       refY="0"
77
       refX="0"
78
       id="Arrow1LendI"
79
       style="overflow:visible">
80
      <path
81
         id="path5515"
82
         d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
83
         style="fill:#808080;fill-rule:evenodd;stroke:#808080;stroke-width:1pt"
84
         transform="matrix(-0.8,0,0,-0.8,-10,0)"
85
         inkscape:connector-curvature="0" />
86
    </marker>
87
    <marker
88
       inkscape:stockid="Arrow1LendI"
89
       orient="auto"
90
       refY="0"
91
       refX="0"
92
       id="Arrow1LendI-2"
93
       style="overflow:visible">
94
      <path
95
         inkscape:connector-curvature="0"
96
         id="path5515-0"
97
         d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
98
         style="fill:#808080;fill-rule:evenodd;stroke:#808080;stroke-width:1pt"
99
         transform="matrix(-0.8,0,0,-0.8,-10,0)" />
100
    </marker>
101
  </defs>
102
  <metadata
103
     id="metadata7">
104
    <rdf:RDF>
105
      <cc:Work
106
         rdf:about="">
107
        <dc:format>image/svg+xml</dc:format>
108
        <dc:type
109
           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
110
        <dc:title></dc:title>
111
      </cc:Work>
112
    </rdf:RDF>
113
  </metadata>
114
  <g
115
     inkscape:label="Layer 1"
116
     inkscape:groupmode="layer"
117
     id="layer1"
118
     transform="translate(-36.278871,-279.59891)">
119
    <path
120
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
121
       d="M 587.04407,408.86412 556.35097,447.755"
122
       id="path6029"
123
       inkscape:connector-type="polyline"
124
       inkscape:connector-curvature="0"
125
       inkscape:connection-start="#rect4439-36"
126
       inkscape:connection-start-point="d4" />
127
    <path
128
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
129
       d="M 519.33173,410.88441 556.35097,447.755"
130
       id="path6031"
131
       inkscape:connector-type="polyline"
132
       inkscape:connector-curvature="0"
133
       inkscape:connection-start="#rect4439-3"
134
       inkscape:connection-start-point="d4" />
135
    <path
136
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
137
       d="m 556.35097,157.4873 0,35.48559"
138
       id="path6033"
139
       inkscape:connector-type="polyline"
140
       inkscape:connector-curvature="0"
141
       inkscape:connection-end="#path4443-3"
142
       inkscape:connection-end-point="d4"
143
       transform="translate(0,308.2677)" />
144
    <path
145
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
146
       d="m 228.12949,609.89385 -106,74.25338"
147
       id="path5754"
148
       inkscape:connector-type="polyline"
149
       inkscape:connector-curvature="0" />
150
    <path
151
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
152
       d="m 227.11934,769.51228 159.92228,74.75846"
153
       id="path4548"
154
       inkscape:connector-type="polyline"
155
       inkscape:connector-curvature="0" />
156
    <path
157
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
158
       d="m 122.12949,684.14723 104.98985,85.36505"
159
       id="path5752"
160
       inkscape:connector-type="polyline"
161
       inkscape:connector-curvature="0" />
162
    <path
163
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
164
       d="m 228.12949,609.89385 -1.01015,159.61843"
165
       id="path5756"
166
       inkscape:connector-type="polyline"
167
       inkscape:connector-curvature="0" />
168
    <path
169
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
170
       d="m 387.04162,844.27074 0,69.70054"
171
       id="path5758"
172
       inkscape:connector-type="polyline"
173
       inkscape:connector-curvature="0" />
174
    <path
175
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
176
       d="m 387.04162,913.97128 0,71.72083"
177
       id="path5760"
178
       inkscape:connector-type="polyline"
179
       inkscape:connector-curvature="0" />
180
    <path
181
       sodipodi:type="arc"
182
       style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0"
183
       id="path4447-4"
184
       sodipodi:cx="165.66502"
185
       sodipodi:cy="351.14514"
186
       sodipodi:rx="91.923882"
187
       sodipodi:ry="26.263966"
188
       d="m 257.58891,351.14514 c 0,14.50519 -41.15573,26.26397 -91.92389,26.26397 -50.76815,0 -91.923878,-11.75878 -91.923878,-26.26397 0,-14.50519 41.155728,-26.26396 91.923878,-26.26396 50.76816,0 91.92389,11.75877 91.92389,26.26396 z"
189
       transform="translate(221.3766,562.82614)" />
190
    <path
191
       sodipodi:type="arc"
192
       style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0"
193
       id="path4447-1"
194
       sodipodi:cx="165.66502"
195
       sodipodi:cy="351.14514"
196
       sodipodi:rx="91.923882"
197
       sodipodi:ry="26.263966"
198
       d="m 257.58891,351.14514 c 0,14.50519 -41.15573,26.26397 -91.92389,26.26397 -50.76815,0 -91.923878,-11.75878 -91.923878,-26.26397 0,-14.50519 41.155728,-26.26396 91.923878,-26.26396 50.76816,0 91.92389,11.75877 91.92389,26.26396 z"
199
       transform="translate(221.3766,634.54697)" />
200
    <path
201
       sodipodi:type="arc"
202
       style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0"
203
       id="path4447-7"
204
       sodipodi:cx="165.66502"
205
       sodipodi:cy="351.14514"
206
       sodipodi:rx="91.923882"
207
       sodipodi:ry="26.263966"
208
       d="m 257.58891,351.14514 c 0,14.50519 -41.15573,26.26397 -91.92389,26.26397 -50.76815,0 -91.923878,-11.75878 -91.923878,-26.26397 0,-14.50519 41.155728,-26.26396 91.923878,-26.26396 50.76816,0 91.92389,11.75877 91.92389,26.26396 z"
209
       transform="translate(221.3766,493.1256)" />
210
    <path
211
       sodipodi:type="arc"
212
       style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0"
213
       id="path4447"
214
       sodipodi:cx="165.66502"
215
       sodipodi:cy="351.14514"
216
       sodipodi:rx="91.923882"
217
       sodipodi:ry="26.263966"
218
       d="m 257.58891,351.14514 c 0,14.50519 -41.15573,26.26397 -91.92389,26.26397 -50.76815,0 -91.923878,-11.75878 -91.923878,-26.26397 0,-14.50519 41.155728,-26.26396 91.923878,-26.26396 50.76816,0 91.92389,11.75877 91.92389,26.26396 z"
219
       transform="translate(61.454315,418.36714)" />
220
    <path
221
       sodipodi:type="arc"
222
       style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-dasharray:1, 1;stroke-dashoffset:0"
223
       id="path4445"
224
       sodipodi:cx="153.03812"
225
       sodipodi:cy="292.05121"
226
       sodipodi:rx="85.357887"
227
       sodipodi:ry="27.779196"
228
       d="m 238.396,292.05121 c 0,15.34202 -38.21602,27.77919 -85.35788,27.77919 -47.14186,0 -85.357891,-12.43717 -85.357891,-27.77919 0,-15.34203 38.216031,-27.7792 85.357891,-27.7792 47.14186,0 85.35788,12.43717 85.35788,27.7792 z"
229
       transform="translate(-30.908629,392.09602)" />
230
    <path
231
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
232
       d="m 275.74209,228.41161 -47.6126,73.21454"
233
       id="path6039"
234
       inkscape:connector-type="polyline"
235
       inkscape:connector-curvature="0"
236
       transform="translate(0,308.2677)" />
237
    <path
238
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
239
       d="m 175.09648,228.4116 53.03301,73.21455"
240
       id="path6037"
241
       inkscape:connector-type="polyline"
242
       inkscape:connector-curvature="0"
243
       transform="translate(0,308.2677)" />
244
    <path
245
       sodipodi:type="arc"
246
       style="fill:#ffffff;stroke:#000000"
247
       id="path4443"
248
       sodipodi:cx="168.69548"
249
       sodipodi:cy="245.07913"
250
       sodipodi:rx="66.670067"
251
       sodipodi:ry="15.152288"
252
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
253
       transform="translate(59.434007,364.81472)" />
254
    <path
255
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
256
       d="m 281.59076,395.18297 -5.84867,141.49634"
257
       id="path5768"
258
       inkscape:connector-type="polyline"
259
       inkscape:connector-curvature="0" />
260
    <path
261
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
262
       d="m 142.89347,465.75497 132.84862,70.92434"
263
       id="path5764"
264
       inkscape:connector-type="polyline"
265
       inkscape:connector-curvature="0" />
266
    <rect
267
       style="fill:#ffffff;stroke:#000000"
268
       id="rect4439-1"
269
       width="72.73098"
270
       height="39.39595"
271
       x="239.3766"
272
       y="516.98132" />
273
    <path
274
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
275
       d="M 281.59076,395.18297 175.09647,536.67931"
276
       id="path5766"
277
       inkscape:connector-type="polyline"
278
       inkscape:connector-curvature="0" />
279
    <path
280
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
281
       d="m 140.69037,395.79245 2.2031,69.96252"
282
       id="path5770"
283
       inkscape:connector-type="polyline"
284
       inkscape:connector-curvature="0" />
285
    <path
286
       sodipodi:type="arc"
287
       style="fill:#ffffff;stroke:#000000;stroke-width:0.61840045;stroke-miterlimit:4;stroke-dasharray:none"
288
       id="path4443-3-8-4"
289
       sodipodi:cx="168.69548"
290
       sodipodi:cy="245.07913"
291
       sodipodi:rx="66.670067"
292
       sodipodi:ry="15.152288"
293
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
294
       transform="matrix(1.3139941,0,0,1.9900637,334.6861,290.16903)" />
295
    <path
296
       sodipodi:type="arc"
297
       style="fill:#ffffff;stroke:#000000;stroke-width:0.50245172;stroke-miterlimit:4;stroke-dasharray:none"
298
       id="path4443-3-8"
299
       sodipodi:cx="168.69548"
300
       sodipodi:cy="245.07913"
301
       sodipodi:rx="66.670067"
302
       sodipodi:ry="15.152288"
303
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
304
       transform="matrix(1.697119,0,0,2.3339904,270.05466,27.535178)" />
305
    <path
306
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
307
       d="m 142.89347,465.75497 32.203,70.92434"
308
       id="path5762"
309
       inkscape:connector-type="polyline"
310
       inkscape:connector-curvature="0" />
311
    <rect
312
       style="fill:none"
313
       id="rect3010"
314
       width="434.28571"
315
       height="74.285713"
316
       x="87.428574"
317
       y="279.59891" />
318
    <text
319
       xml:space="preserve"
320
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
321
       x="156.27484"
322
       y="349.56943"
323
       id="text3782"
324
       sodipodi:linespacing="125%"><tspan
325
         style="font-size:16px"
326
         sodipodi:role="line"
327
         id="tspan3784"
328
         x="156.27484"
329
         y="349.56943">MODIS data</tspan></text>
330
    <text
331
       xml:space="preserve"
332
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
333
       x="321.29907"
334
       y="496.00278"
335
       id="text2987-0-8"
336
       sodipodi:linespacing="125%"><tspan
337
         style="font-size:12px;text-align:center;text-anchor:middle"
338
         id="tspan3823-7"
339
         sodipodi:role="line"
340
         x="321.29907"
341
         y="496.00278" /></text>
342
    <text
343
       xml:space="preserve"
344
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
345
       x="546.98883"
346
       y="348.46786"
347
       id="text2987-4"
348
       sodipodi:linespacing="125%"><tspan
349
         style="font-size:16px;text-align:center;text-anchor:middle"
350
         sodipodi:role="line"
351
         x="546.98883"
352
         y="348.46786"
353
         id="tspan4124">Meteorological Station Data</tspan></text>
354
    <text
355
       xml:space="preserve"
356
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
357
       x="570.02509"
358
       y="674.41626"
359
       id="text2987-4-53-9-5-0"
360
       sodipodi:linespacing="125%"><tspan
361
         style="font-size:12px;text-align:center;text-anchor:middle"
362
         sodipodi:role="line"
363
         x="570.02509"
364
         y="674.41626"
365
         id="tspan4271-1" /></text>
366
    <text
367
       xml:space="preserve"
368
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
369
       x="330.39041"
370
       y="837.4798"
371
       id="text2987-4-53-9-5-9-2"
372
       sodipodi:linespacing="125%"><tspan
373
         style="font-size:12px;text-align:center;text-anchor:middle"
374
         sodipodi:role="line"
375
         x="330.39041"
376
         y="837.4798"
377
         id="tspan4271-2-3" /></text>
378
    <g
379
       id="g6130">
380
      <rect
381
         y="372.55896"
382
         x="84.121834"
383
         height="46.467018"
384
         width="113.13708"
385
         id="rect4435"
386
         style="fill:#ffffff;stroke:#000000" />
387
      <text
388
         sodipodi:linespacing="125%"
389
         id="text2987"
390
         y="393.14954"
391
         x="140.49994"
392
         style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
393
         xml:space="preserve"><tspan
394
           y="393.14954"
395
           x="140.49994"
396
           id="tspan2989"
397
           sodipodi:role="line"
398
           style="font-size:14px;text-align:center;text-anchor:middle">MOD06_L2</tspan><tspan
399
           id="tspan3835"
400
           y="408.65604"
401
           x="140.49994"
402
           sodipodi:role="line"
403
           style="font-size:12px;text-align:center;text-anchor:middle">Cloud Product</tspan></text>
404
    </g>
405
    <g
406
       id="g6124">
407
      <rect
408
         y="371.23422"
409
         x="221.21829"
410
         height="47.897472"
411
         width="120.74493"
412
         id="rect4437"
413
         style="fill:#ffffff;stroke:#000000;stroke-width:0.78849989" />
414
      <text
415
         sodipodi:linespacing="125%"
416
         id="text2987-0"
417
         y="391.37698"
418
         x="280.96957"
419
         style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
420
         xml:space="preserve"><tspan
421
           y="391.37698"
422
           x="280.96957"
423
           id="tspan2989-4"
424
           sodipodi:role="line"
425
           style="font-size:14px">MOD11A1</tspan><tspan
426
           y="406.88348"
427
           x="280.96957"
428
           sodipodi:role="line"
429
           id="tspan3823"
430
           style="font-size:12px;text-align:center;text-anchor:middle">Land Surface Temp</tspan></text>
431
    </g>
432
    <rect
433
       style="fill:#ffffff;stroke:#000000"
434
       id="rect4439"
435
       width="72.73098"
436
       height="39.39595"
437
       x="138.73099"
438
       y="516.98132" />
439
    <text
440
       xml:space="preserve"
441
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
442
       x="175.05545"
443
       y="524.98987"
444
       id="text2987-0-4-8-3"
445
       sodipodi:linespacing="125%"><tspan
446
         style="font-size:12px;text-align:center;text-anchor:middle"
447
         sodipodi:role="line"
448
         x="175.05545"
449
         y="524.98987"
450
         id="tspan3912-5" /><tspan
451
         style="font-size:12px;text-align:center;text-anchor:middle"
452
         sodipodi:role="line"
453
         x="175.05545"
454
         y="539.98987"
455
         id="tspan3946">QA Flags</tspan></text>
456
    <g
457
       id="g6118">
458
      <text
459
         sodipodi:linespacing="125%"
460
         id="text2987-4-9"
461
         y="388.05463"
462
         x="499.36975"
463
         style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
464
         xml:space="preserve"><tspan
465
           id="tspan4124-6"
466
           y="388.05463"
467
           x="499.36975"
468
           sodipodi:role="line"
469
           style="font-size:12px;text-align:center;text-anchor:middle">GHCND</tspan><tspan
470
           id="tspan4152"
471
           y="403.05463"
472
           x="499.36975"
473
           sodipodi:role="line"
474
           style="font-size:12px;text-align:center;text-anchor:middle">dataset</tspan></text>
475
      <rect
476
         y="371.48846"
477
         x="463.18884"
478
         height="39.39595"
479
         width="72.73098"
480
         id="rect4439-3"
481
         style="fill:none;stroke:#000000" />
482
    </g>
483
    <g
484
       id="g6112">
485
      <text
486
         sodipodi:linespacing="125%"
487
         id="text2987-4-5"
488
         y="386.1398"
489
         x="602.55475"
490
         style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
491
         xml:space="preserve"><tspan
492
           id="tspan4167"
493
           y="386.1398"
494
           x="602.55475"
495
           sodipodi:role="line"
496
           style="font-size:12px;text-align:center;text-anchor:middle">Other</tspan><tspan
497
           id="tspan4171"
498
           y="401.1398"
499
           x="602.55475"
500
           sodipodi:role="line"
501
           style="font-size:12px;text-align:center;text-anchor:middle">databases</tspan></text>
502
      <rect
503
         y="369.46817"
504
         x="566.22443"
505
         height="39.39595"
506
         width="72.73098"
507
         id="rect4439-36"
508
         style="fill:none;stroke:#000000" />
509
    </g>
510
    <g
511
       id="g6136"
512
       transform="translate(0,8)">
513
      <rect
514
         style="fill:#ffffff;stroke:#000000;stroke-width:1;stroke-miterlimit:3.9000001;stroke-opacity:1;stroke-dasharray:none;stroke-dashoffset:0"
515
         id="rect5799"
516
         width="151.52289"
517
         height="34.345188"
518
         x="480.58951"
519
         y="430.5824" />
520
      <text
521
         sodipodi:linespacing="125%"
522
         id="text2987-4-3"
523
         y="451.06555"
524
         x="556.08142"
525
         style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
526
         xml:space="preserve"><tspan
527
           id="tspan4124-9"
528
           y="451.06555"
529
           x="556.08142"
530
           sodipodi:role="line"
531
           style="font-size:12px;text-align:center;text-anchor:middle">PostgreSQL database</tspan></text>
532
    </g>
533
    <text
534
       xml:space="preserve"
535
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
536
       x="283.89996"
537
       y="534.44983"
538
       id="text2987-0-4-8-3-7"
539
       sodipodi:linespacing="125%"><tspan
540
         style="font-size:12px;text-align:center;text-anchor:middle"
541
         sodipodi:role="line"
542
         x="283.89996"
543
         y="534.44983"
544
         id="tspan3974" /></text>
545
    <text
546
       xml:space="preserve"
547
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
548
       x="228.33749"
549
       y="613.20441"
550
       id="text2987-0-4-8-3-3-4-7"
551
       sodipodi:linespacing="125%"><tspan
552
         style="font-size:12px;text-align:center;text-anchor:middle"
553
         sodipodi:role="line"
554
         x="228.33749"
555
         y="613.20441"
556
         id="tspan4017-4">QA Screening</tspan></text>
557
    <text
558
       xml:space="preserve"
559
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
560
       x="121.97128"
561
       y="672.45776"
562
       id="text2987-0-4-8-3-3-4"
563
       sodipodi:linespacing="125%"><tspan
564
         style="font-size:12px;text-align:center;text-anchor:middle;fill:#4d4d4d"
565
         sodipodi:role="line"
566
         x="121.97128"
567
         y="672.45776"
568
         id="tspan3989-7">MOD06 only</tspan><tspan
569
         style="font-size:12px;text-align:center;text-anchor:middle;fill:#4d4d4d"
570
         sodipodi:role="line"
571
         x="121.97128"
572
         y="687.45776"
573
         id="tspan4096">Generate daily</tspan><tspan
574
         style="font-size:12px;text-align:center;text-anchor:middle;fill:#4d4d4d"
575
         sodipodi:role="line"
576
         x="121.97128"
577
         y="702.45776"
578
         id="tspan4017">product</tspan></text>
579
    <text
580
       xml:space="preserve"
581
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
582
       x="226.96114"
583
       y="765.32288"
584
       id="text2987-0-4-8-1"
585
       sodipodi:linespacing="125%"><tspan
586
         style="font-size:12px;text-align:center;text-anchor:middle"
587
         sodipodi:role="line"
588
         x="226.96114"
589
         y="765.32288"
590
         id="tspan3912-2">Calculate monthly</tspan><tspan
591
         style="font-size:12px;text-align:center;text-anchor:middle"
592
         sodipodi:role="line"
593
         x="226.96114"
594
         y="780.32288"
595
         id="tspan4055">climatologies</tspan></text>
596
    <text
597
       xml:space="preserve"
598
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
599
       x="387.41956"
600
       y="832.5813"
601
       id="text2987-0-4-8"
602
       sodipodi:linespacing="125%"><tspan
603
         style="font-size:12px;text-align:center;text-anchor:middle"
604
         sodipodi:role="line"
605
         x="387.41956"
606
         y="832.5813"
607
         id="tspan3863-9">Reproject to</tspan><tspan
608
         style="font-size:12px;text-align:center;text-anchor:middle"
609
         sodipodi:role="line"
610
         x="389.32971"
611
         y="847.5813"
612
         id="tspan3910">local projection </tspan><tspan
613
         style="font-size:12px;text-align:center;text-anchor:middle"
614
         sodipodi:role="line"
615
         x="387.41956"
616
         y="862.5813"
617
         id="tspan3912">for interpolation</tspan></text>
618
    <text
619
       xml:space="preserve"
620
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
621
       x="387.04453"
622
       y="909.78186"
623
       id="text2987-4-53-9-5-9"
624
       sodipodi:linespacing="125%"><tspan
625
         style="font-size:12px;text-align:center;text-anchor:middle"
626
         sodipodi:role="line"
627
         x="387.04453"
628
         y="909.78186"
629
         id="tspan4348">Daily interpolation</tspan><tspan
630
         style="font-size:12px;text-align:center;text-anchor:middle"
631
         sodipodi:role="line"
632
         x="387.04453"
633
         y="924.78186"
634
         id="tspan4346">(Kriging, GAM, GWR)</tspan></text>
635
    <text
636
       xml:space="preserve"
637
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
638
       x="387.04456"
639
       y="981.50269"
640
       id="text2987-4-53-9-5-9-1"
641
       sodipodi:linespacing="125%"><tspan
642
         style="font-size:12px;text-align:center;text-anchor:middle"
643
         sodipodi:role="line"
644
         x="387.04456"
645
         y="981.50269"
646
         id="tspan4348-8">Daily Predictions</tspan><tspan
647
         style="font-size:12px;text-align:center;text-anchor:middle"
648
         sodipodi:role="line"
649
         x="387.04456"
650
         y="996.50269"
651
         id="tspan4346-9">(Kriging, GAM, GWR)</tspan></text>
652
    <path
653
       sodipodi:type="arc"
654
       style="fill:#ffffff;stroke:#000000"
655
       id="path4443-8"
656
       sodipodi:cx="168.69548"
657
       sodipodi:cy="245.07913"
658
       sodipodi:rx="66.670067"
659
       sodipodi:ry="15.152288"
660
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
661
       transform="matrix(1,0,0,1.3831658,-25.80201,126.7699)" />
662
    <g
663
       id="g4520"
664
       transform="translate(74.70817,-5.26443)">
665
      <text
666
         sodipodi:linespacing="125%"
667
         id="text2987-0-4-8-3-3"
668
         y="466.82996"
669
         x="68.364006"
670
         style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
671
         xml:space="preserve"><tspan
672
           id="tspan3946-6"
673
           y="466.82996"
674
           x="68.364006"
675
           sodipodi:role="line"
676
           style="font-size:12px;text-align:center;text-anchor:middle">Convert swath</tspan><tspan
677
           id="tspan3989"
678
           y="481.82996"
679
           x="68.364006"
680
           sodipodi:role="line"
681
           style="font-size:12px;text-align:center;text-anchor:middle">to grid</tspan></text>
682
    </g>
683
    <text
684
       xml:space="preserve"
685
       style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
686
       x="275.44913"
687
       y="533.46838"
688
       id="text2987-0-4-8-3-9"
689
       sodipodi:linespacing="125%"><tspan
690
         style="font-size:12px;text-align:center;text-anchor:middle"
691
         sodipodi:role="line"
692
         x="275.44913"
693
         y="533.46838"
694
         id="tspan3946-4">Parameters</tspan><tspan
695
         id="tspan5040"
696
         style="font-size:12px;text-align:center;text-anchor:middle"
697
         sodipodi:role="line"
698
         x="275.44913"
699
         y="548.46838">of Interest</tspan></text>
700
    <path
701
       sodipodi:type="arc"
702
       style="fill:#ffffff;stroke:#000000"
703
       id="path4443-3"
704
       sodipodi:cx="168.69548"
705
       sodipodi:cy="245.07913"
706
       sodipodi:rx="66.670067"
707
       sodipodi:ry="15.152288"
708
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
709
       transform="matrix(1.1419854,0,0,1.3914026,363.7032,181.31978)" />
710
    <path
711
       sodipodi:type="arc"
712
       style="fill:#ffffff;stroke:#000000;stroke-width:0.50245172;stroke-miterlimit:4;stroke-dasharray:none"
713
       id="path4443-3-8-0"
714
       sodipodi:cx="168.69548"
715
       sodipodi:cy="245.07913"
716
       sodipodi:rx="66.670067"
717
       sodipodi:ry="15.152288"
718
       d="m 235.36555,245.07913 c 0,8.36838 -29.84921,15.15229 -66.67007,15.15229 -36.82086,0 -66.67007,-6.78391 -66.67007,-15.15229 0,-8.36838 29.84921,-15.15229 66.67007,-15.15229 36.82086,0 66.67007,6.78391 66.67007,15.15229 z"
719
       transform="matrix(1.697119,0,0,2.3339904,270.05466,119.40639)" />
720
    <text
721
       xml:space="preserve"
722
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
723
       x="556.30994"
724
       y="773.7027"
725
       id="text2987-4-53-9-5-9-8"
726
       sodipodi:linespacing="125%"><tspan
727
         style="font-size:12px;text-align:center;text-anchor:middle"
728
         sodipodi:role="line"
729
         x="556.30994"
730
         y="773.7027"
731
         id="tspan4346-8">Extraction of data</tspan><tspan
732
         style="font-size:12px;text-align:center;text-anchor:middle"
733
         sodipodi:role="line"
734
         x="556.30994"
735
         y="788.7027"
736
         id="tspan4380">for interpolation</tspan></text>
737
    <text
738
       xml:space="preserve"
739
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
740
       x="556.57361"
741
       y="518.13409"
742
       id="text2987-4-53"
743
       sodipodi:linespacing="125%"><tspan
744
         style="font-size:12px;text-align:center;text-anchor:middle"
745
         sodipodi:role="line"
746
         x="556.57361"
747
         y="518.13409"
748
         id="tspan4124-4">Apply</tspan><tspan
749
         style="font-size:12px;text-align:center;text-anchor:middle"
750
         sodipodi:role="line"
751
         x="556.57361"
752
         y="533.13409"
753
         id="tspan4199">Quality Control</tspan></text>
754
    <text
755
       xml:space="preserve"
756
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
757
       x="556.35388"
758
       y="583.85809"
759
       id="text2987-4-53-9"
760
       sodipodi:linespacing="125%"><tspan
761
         style="font-size:12px;text-align:center;text-anchor:middle"
762
         sodipodi:role="line"
763
         x="558.26404"
764
         y="583.85809"
765
         id="tspan4199-9">Annotate station </tspan><tspan
766
         style="font-size:12px;text-align:center;text-anchor:middle"
767
         sodipodi:role="line"
768
         x="556.35388"
769
         y="598.85809"
770
         id="tspan4227">data with covariates</tspan><tspan
771
         style="font-size:12px;text-align:center;text-anchor:middle"
772
         sodipodi:role="line"
773
         x="556.35388"
774
         y="613.85809"
775
         id="tspan4229">(elevation, slope, aspect, etc)</tspan></text>
776
    <text
777
       xml:space="preserve"
778
       style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
779
       x="556.30994"
780
       y="687.22931"
781
       id="text2987-4-53-9-5"
782
       sodipodi:linespacing="125%"><tspan
783
         style="font-size:12px;text-align:center;text-anchor:middle"
784
         sodipodi:role="line"
785
         x="556.30994"
786
         y="687.22931"
787
         id="tspan4267">Generate monthly climatologies</tspan><tspan
788
         style="font-size:12px;text-align:center;text-anchor:middle"
789
         sodipodi:role="line"
790
         x="556.30994"
791
         y="702.22931"
792
         id="tspan4271">for 'climate-aided' interpolation</tspan></text>
793
    <path
794
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
795
       d="m 556.35097,418.51633 0,20.95407"
796
       id="path6017"
797
       inkscape:connector-type="orthogonal"
798
       inkscape:connector-curvature="0"
799
       inkscape:connection-start="#path4443-3-8-0"
800
       inkscape:connection-start-point="d4"
801
       inkscape:connection-end="#path4443-3-8-4"
802
       inkscape:connection-end-point="d4"
803
       transform="translate(0,308.2677)" />
804
    <path
805
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
806
       d="m 556.35098,235.13876 0,20.77576"
807
       id="path6019"
808
       inkscape:connector-type="orthogonal"
809
       inkscape:connector-curvature="0"
810
       inkscape:connection-start="#path4443-3"
811
       inkscape:connection-start-point="d4"
812
       inkscape:connection-end="#path4443-3-8"
813
       inkscape:connection-end-point="d4"
814
       transform="translate(0,308.2677)" />
815
    <path
816
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
817
       d="m 556.35096,326.64512 0,21.14062"
818
       id="path6021"
819
       inkscape:connector-type="orthogonal"
820
       inkscape:connector-curvature="0"
821
       inkscape:connection-start="#path4443-3-8"
822
       inkscape:connection-start-point="d4"
823
       inkscape:connection-end="#path4443-3-8-0"
824
       inkscape:connection-end-point="d4"
825
       transform="translate(0,308.2677)" />
826
    <path
827
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
828
       d="m 498.55233,492.28467 -57.36768,22.49131"
829
       id="path6110"
830
       inkscape:connector-type="polyline"
831
       inkscape:connector-curvature="0"
832
       inkscape:connection-start="#path4443-3-8-4"
833
       inkscape:connection-start-point="d4"
834
       inkscape:connection-end="#path4447-7"
835
       inkscape:connection-end-point="d4"
836
       transform="translate(0,308.2677)" />
837
  </g>
838
</svg>
climate/extra/landcover_lst.R
1
### Script to explore the relationships between LST and land cover in Oregon
2
### The goal is to inform the interpolation of LST using land cover
3

  
4

  
climate/procedures/MOD06_L2_data_summary.r
39 39
##########################
40 40
#### explore the data
41 41

  
42
## load data
43 42
months=seq(as.Date("2000-01-15"),as.Date("2000-12-15"),by="month")
43

  
44

  
45
## load data
44 46
cerfiles=list.files(summarydatadir,pattern="CER_mean_.*tif$",full=T); cerfiles
45 47
cer=brick(stack(cerfiles))
46 48
setZ(cer,months,name="time")
......
71 73
cldm=mean(cld,na.rm=T)
72 74
### TODO: change to bilinear if reprojecting!
73 75

  
76
cer20files=list.files(summarydatadir,pattern="CER_P20um_.*tif$",full=T); cer20files
77
cer20=brick(stack(cer20files))
78
setZ(cer20,months,name="time")
79
cer20@z=list(months)
80
cer20@zname="time"
81
layerNames(cer20) <- as.character(format(months,"%b"))
82
#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")
83
cotm=mean(cot,na.rm=T)
84
### TODO: change to bilinear!
85

  
86

  
74 87
### load PRISM data for comparison
75 88
prism=brick("data/prism/prism_climate.nc",varname="ppt")
76 89
## project to sinusoidal
......
83 96
layerNames(prism) <- as.character(format(months,"%b"))
84 97

  
85 98
####  build a pixel by variable matrix
86
vars=c("cer","cld","cot","prism")
99
vars=c("cer","cer20","cld","cot","prism")
87 100
bd=melt(as.matrix(vars[1]))
88 101
colnames(bd)=c("cell","month",vars[1])
89 102
for(v in vars[-1]) {print(v); bd[,v]=melt(as.matrix(get(v)))$value}
......
116 129

  
117 130
### extract MOD06 data for each station
118 131
stcer=extract(cer,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
132
stcer20=extract(cer20,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
119 133
stcot=extract(cot,st2)#;colnames(stcot)=paste("cot_mean_",1:12,sep="")
120 134
stcld=extract(cld,st2)#;colnames(stcld)=paste("cld_mean_",1:12,sep="")
121
mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcot,stcld)
135
mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcer20,stcot,stcld)
122 136
mod06l=melt(mod06,id.vars=c("id","lon","lat"))
123 137
mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
124 138
mod06l=as.data.frame(cast(mod06l,id+lon+lat+month~variable+moment,value="value"))
......
158 172

  
159 173
## melt it
160 174
mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))
161
levels(mod06sl$variable)=c("Effective Radius (um)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
175
levels(mod06sl$variable)=c("Effective Radius (um)","Very Cloudy Days (%)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
162 176

  
163 177
###################################################################
164 178
###################################################################
......
234 248
print(p)
235 249

  
236 250
### monthly comparisons of variables
237
mod06sl=melt(mod06s,measure.vars=c("value","COT_mean","CER_mean"))
238
bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
239
splom(mod06s[grep("CER|COT|CLD",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
251
#mod06sl=melt(mod06s,measure.vars=c("ppt","COT_mean","CER_mean","CER_P20um"))
252
#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
253
#splom(mod06s[grep("CER|COT|CLD|ppt",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
240 254

  
241 255
### run some regressions
242 256
#plot(log(ppt)~COT_mean,data=mod06s)
......
246 260
 xyplot(ppt~value|variable,groups=glon,data=mod06sl,
247 261
       scales=list(y=list(log=T),x=list(relation="free",log=F)),
248 262
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"),
249
       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(4,1))+
263
       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(5,1))+
250 264
  layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
251
  layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
265
  layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)+
266
  layer(panel.abline(lm(y~x),col="red"))+
267
  layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
268

  
252 269

  
253 270
## ppt~metric with longitude bins
254 271
#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)
......
257 274
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
258 275
p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
259 276
   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black'))
277
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
278

  
260 279
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
261 280
p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
262 281
    )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black'))
......
268 287
  layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
269 288
  layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
270 289

  
271
save(p1,p2,p3,file="plotdata.Rdata")
272
load("plotdata.Rdata")
290
p4=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Very Cloudy Days (%)",],
291
       scales=list(y=list(log=T),x=list(relation="free",log=F)),
292
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
293
ylab="Mean Monthly Station Precipitation (mm)",xlab="Proportion days with Cloud Effective Radius >20um from MOD06_L2 (%)",layout=c(1,1))+
294
  layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
295
  layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
296

  
297
#save(p1,p2,p3,file="plotdata.Rdata")
298
#load("plotdata.Rdata")
273 299

  
274
CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
300
#CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
275 301
print(p3,position=c(0,0,1,.5),save.object=F)
302
#print(p4,position=c(0,0,1,.5),save.object=F)
276 303
print(p1,split=c(1,1,2,2),new=F)
277 304
print(p2,split=c(2,1,2,2),new=F)
278
dev.off()
279
system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
305
#dev.off()
306
#system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
280 307
                                        #dev.off()
281 308

  
282 309
## with elevation
climate/procedures/interpolation.r
1
#! /bin/r
2

  
3

  
4
#system("ssh -X turaco")
5
#export LD_LIBRARY_PATH=${LD_LIBRARY_PATH:+$LD_LIBRARY_PATH:}/home/adamw/local/lib
6
#export KMP_DUPLICATE_LIB_OK=TRUE
7
#R
8
### Download and clean up GHCN data for CFR
9
#install.packages("ncdf4",lib="/home/adamw/rlib",configure.args="--with-nc-config=/home/adamw/local/bin/nc-config")
10

  
11
library(sp);library(rgdal);
12
library(reshape)
13
library(ncdf4)
14
library(geosphere)
15
library(rgeos)
16
library(multicore);library(RSQLite)
17
library(spBayes)
18
library(geoR) #for conventional kriging
19
library(raster)
20
## for review/analysis
21
library(coda)
22
library(lattice)
23
library(MBA)
24

  
25
###  Load functions
26
#source("code/interpolationfunc.r")
27

  
28
X11.options(type="Xlib")
29
ncores=20  #number of threads to use
30

  
31
setwd("/home/adamw/acrobates/projects/interp")
32

  
33
### read in full dataset
34
load("data/ghcn/roi_ghcn.Rdata")
35
load("stroi.Rdata")
36
source("code/GHCN_functions.r")
37

  
38
## rescale units to C and mm
39
d$value=d$value/10.0
40

  
41
## add flag for Benoit's training stations to identify training/validation
42
stt=readOGR("data/ghcn/tempdata/ghcn_1507_s.shp","ghcn_1507_s")$STAT_ID
43
d$fit=d$id%in%stt
44

  
45
### read in and subset to Benoit's dates
46
dates=as.Date(as.character(
47
  read.table("data/ghcn/tempdata/dates_interpolation_03052012.txt")$V1),"%Y%m%d")
48

  
49
##########################################################
50
#### load region of interest
51
  roi=readOGR("data/boundaries/statesp020.shp","statesp020")
52
  proj4string(roi)=CRS("+proj=longlat +datum=WGS84")
53
  roi=roi[roi$STATE=="Oregon",]
54
  roi=gUnionCascaded(roi)  #dissolve any subparts of roi
55
  
56
  ## distance from ROI to buffer for border stations
57
  dis=100 #km
58
  ## buffer roi (transform to azimuthal equidistant with centroid of roi, add 'dis' buffer, then transform back)
59
  roic=centroid(roi)
60
  
61
  roib=spTransform(
62
    gBuffer(spTransform(roi,
63
                        CRS(paste("+proj=aeqd +lat_0=",roic[2]," +lon_0=",
64
                                  roic[1]," +ellps=WGS84 +datum=WGS84 +units=m +no_defs ",sep=""))),
65
            width=dis*1000),
66
    CRS("+proj=longlat +datum=WGS84"))
67

  
68
### generate knots
69
knots=makegrid(roib,cellsize=0.5)
70
coordinates(knots)=c("x1","x2")
71
proj4string(knots)=proj4string(roi)
72
knots=knots[!is.na(overlay(knots,roi))]
73

  
74

  
75
###########################################################3
76
#### interpolation function
77

  
78
interp<-function(parms,niter=100,istart=60,ithin=1){
79
  ## load variables from parameter list
80
  idate=as.Date(parms$date)
81
  ivar=parms$variable
82
  type=parms$type
83
  model=parms$model
84
  
85
  ## subset single date from data
86
  td=d[d$date==idate,]
87
  td=data.frame(cast(td,id+fit~variable,value="value"))
88
  td$elev=stroi$elev[match(td$id,stroi$id)]
89
  td[c("lon","lat")]=stroi@data[match(td$id,stroi$id),c("lon","lat")]
90
  td$intercept=1
91

  
92
  ## drop any points without coordinates
93
  td=td[!is.na(td$lat)&!is.na(td$lon),]
94
  
95
###########################
96
  ## apply climate correction
97
  if(type=="anom"){
98
    cppt=brick("data/prism/prism_climate.nc",varname="ppt")
99
    cppt=subset(cppt,subset=as.numeric(format(idate,"%m")))
100
    cppt[cppt==-99.99]=NA
101
    td$cppt=as.numeric(extract(cppt,td[,c("lon","lat")],method="bilinear"))
102
    ## Tmax
103
    ctmax=brick("data/prism/prism_climate.nc",varname="tmax")
104
    ctmax=subset(ctmax,subset=as.numeric(format(idate,"%m")))
105
    ctmax[ctmax==-99.99]=NA
106
    td$ctmax=as.numeric(extract(ctmax,td[,c("lon","lat")],method="bilinear"))
107
    ## Tmin
108
    ctmin=brick("data/prism/prism_climate.nc",varname="tmin")
109
    ctmin=subset(ctmin,subset=as.numeric(format(idate,"%m")))
110
    ctmin[ctmin==-99.99]=NA
111
    td$ctmin=as.numeric(extract(ctmin,td[,c("lon","lat")],method="bilinear"))
112
    ## Calculate anomalies
113
    td$atmax=td$tmax-td$ctmax
114
    td$atmin=td$tmin-td$ctmin
115
    ## Scale for precipitation (add 1 to climate to avoid dividing by 0, will be subtracted off later)
116
    ## Add 1 to final value to avoid transformation troubles (geoR will only transform if >0
117
    td$appt=(td$ppt/(ifelse(td$cppt<0,0,td$cppt)+1))+1
118
  }
119

  
120
### Add region flag to data
121
  td$region=ifelse(!is.na(overlay(SpatialPoints(td[,c("lon","lat")]),roi)),"roi","roib")
122
  td=td[!is.na(overlay(SpatialPoints(td[,c("lon","lat")]),roib)),]
123

  
124
#### Extract data for validation
125
fd=td[td$fit&td$region=="roi",]  #fitting dataset 
126
pd=td[!td$fit&td$region=="roi",] #training dataset
127
od=td[!td$fit&td$region!="roi",] #stations outside roi
128

  
129

  
130
### drop missing values
131
fd=fd[!is.na(fd[,as.character(ivar)]),]  #&!is.na(fd$tmin)
132

  
133
### define the formula
134
if(type=="raw"&ivar=="tmax")  f1=formula(tmax~intercept)
135
if(type=="anom"&ivar=="tmax"&model=="intercept")  f1=formula(atmax~intercept)
136
if(type=="anom"&ivar=="tmax"&model=="full")  f1=formula(atmax~lon+lat+elev)
137
  
138
### Run the model
139
  t1=Sys.time()
140
  m1=spLM(f1,data=fd,coords=as.matrix(fd[,c("lon","lat")]), 
141
    starting=list("phi"=2.5,"sigma.sq"=4, "tau.sq"=1),
142
    sp.tuning=list("phi"=0.8, "sigma.sq"=0.1, "tau.sq"=0.2),
143
    priors=list("phi.Unif"=c(0.1, 4), "sigma.sq.IG"=c(2, 1),
144
      "tau.sq.IG"=c(2, 1)),
145
    cov.model="exponential",
146
    knots=coordinates(knots),
147
    n.samples=niter, verbose=TRUE, n.report=1000)
148
  t2=Sys.time()
149

  
150
  
151
#m1s=mcmc(m1$p.samples)
152
#summary(m1s)
153
#xyplot(m1s,scales=list(y=list(rot=0)),main="Posterior Samples of Model Parameters")
154

  
155
### Run the model
156
#q=3
157
#nltr=q*(q-1)/2+q
158

  
159
#m1=spMvLM(list(atmin~1,atmax~1),data=fd,coords=as.matrix(fd[,c("lon","lat")]), 
160
#  starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
161
#  sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
162
#  priors=list("phi.Unif"=c(0.1, 3), "sigma.sq.IG"=c(2, 1),
163
#    "tau.sq.IG"=c(2, 1)),
164
#  cov.model="exponential",
165
#  knots=coordinates(knots),
166
#  n.samples=100, verbose=TRUE, n.report=100)
167

  
168
#save to play with later
169
#save(m1,file="m1.Rdata")
170

  
171

  
172
#load("m1.Rdata")
173

  
174
##############################
175
#### make predictions
176

  
177
### Prediction grid
178
#  pgrid=data.frame(coordinates(ctmax))
179
#  ## crop to bbox
180
#  pgrid=pgrid[pgrid$x>=bbox(roi)[1,1]&pgrid$x<=bbox(roi)[1,2]
181
#    &pgrid$y<=bbox(roi)[2,2]&pgrid$y>=bbox(roi)[2,1],]
182
#  ##  Crop to ROI polygon?
183
#  #pgrid=pgrid[!is.na(overlay(SpatialPoints(pgrid),roi)),]
184
#  rownames(pgrid)=1:nrow(pgrid)
185
#  ncluster=200
186
#  pgrid$cluster=cut(as.numeric(rownames(pgrid)),ncluster,labels=1:ncluster)
187
#  ## predict for the full grid
188
#  m1pg=lapply(1:ncluster,function(i){
189
#    print(paste("Starting cluster:",i))
190
#    m1pg=spPredict(m1,pred.coords=pgrid[pgrid$cluster==i,c("x","y")],
191
#              pred.covars=cbind(intercept=rep(1,nrow(pgrid[pgrid$cluster==i,]))),start=istart,thin=ithin)
192
#    ## Generate summaries of y.pred
193
#    ty=t(apply(m1pg$y.pred,1,function(i) c(mean=mean(i),sd=sd(i),quantile(i,c(0.025,0.5,0.975)))))
194
#    colnames(ty)=paste("y.",c("mean","sd","Q2.5","Q50","Q97.5"),sep="")
195
 #   ## Generate summaries of w.pred
196
  #  tw=t(apply(m1pg$w.pred,1,function(i) c(mean=mean(i),sd=sd(i),quantile(i,c(0.025,0.5,0.975)))))
197
  #  colnames(tw)=paste("w.",c("mean","sd","Q2.5","Q50","Q97.5"),sep="")
198
#    return(cbind(ty,tw))
199
#  })
200
#  m1pgs=do.call(rbind.data.frame,m1pg)
201
  
202
  
203
########### Predict only to validation locations
204
m1p=spPredict(m1,pred.coords=pd[,c("lon","lat")],pred.covars=model.matrix(lm(f1,data=pd)),start=istart,thin=ithin)#, 
205
  
206
m1p.y=mcmc(t(m1p$y.pred),start=istart,thin=ithin)
207
m1p.w=mcmc(t(m1$sp.effects.knots),start=istart,thin=ithin)
208
m1p.ys=t(apply(m1p.y,2,quantile,c(0.025,0.5,0.975),na.rm=T))
209
m1p.ws=t(apply(m1p.w,2,quantile,c(0.025,0.5,0.975),na.rm=T))
210

  
211
### recover original scale if necessary
212
 resp=as.character(attr(terms.formula(f1),"variables"))[2]
213
  ## if modeling anomalies, add the climate back on
214
  if(type=="anom") {
215
    pd2=cbind.data.frame(pd,aQ2.5=m1p.ys[,"2.5%"],aQ50=m1p.ys[,"50%"],aQ97.5=m1p.ys[,"97.5%"])
216
# FIXME!  ctmax below needs to be made general for other variables!
217
    pd2$Q2.5=pd2$ctmax+pd2$aQ2.5
218
    pd2$Q50=pd2$ctmax+pd2$aQ50
219
    pd2$Q97.5=pd2$ctmax+pd2$aQ97.5
220
  }
221
  ## if modeling raw values, don't add the climate back on
222
    if(type!="anom") {
223
    pd2=cbind.data.frame(pd,Q2.5=m1p.ys[,"2.5%"],Q50=m1p.ys[,"50%"],Q97.5=m1p.ys[,"97.5%"])
224
  }
225

  
226

  
227
    ### save model output
228
  save(m1,m1p,pd2,file=paste("output/modeloutput_",
229
                gsub("-","",idate),"_",ivar,"_",type,"_",model,".Rdata",sep=""))
230

  
231

  
232

  
233
####################################
234
#### Accuracy Assement
235
m1.dic=spDiag(m1,start=istart,thin=ithin);m1.dic
236
m1.accuracy=accuracy(pd2[,as.character(ivar)],pd2$Q50,ppt=ivar=="ppt")$summary
237

  
238
  #### Write summary file
239
  m1.summary=cbind.data.frame(date=idate,model=model,var=ivar,type=type,
240
    formula=paste(terms(f1),collapse=" "),
241
    t(m1.dic$DIC),t(m1.accuracy),n.iter=niter,runtime.hours=as.numeric(difftime(t1,t2,units="hours")))
242
  write.csv(m1.summary,file=paste("output/modeloutput_summary_",gsub("-","",idate),"_",ivar,"_",type,".csv",sep=""))
243
  
244
##################################################
245
### some summary plots
246
pdf(width=11,height=8.5,file=paste("output/modeloutput_",gsub("-","",idate),"_",ivar,"_",type,".pdf",sep=""))
247

  
248
## Show knots and stations
249
plot(roib);axis(1);axis(2)
250
plot(roi,add=T)
251
#points(knots,pch=3,cex=.5)
252
points(unique(fd[,c('lon','lat')]),cex=.8,pch=16,col="red")
253
  points(unique(pd[,c('lon','lat')]),cex=.5,pch=16,col="blue")
254
points(unique(od[,c('lon','lat')]),cex=.5,pch=16,col="black")
255
  legend(-118,42,legend=c("Fitting","Validation","Border Stations"),pch=c(16,16,16),col=c("red","blue","black"),bg="white")
256
  
257
### Chain behavior
258
  m1s=mcmc(m1$p.samples)
259
  print(xyplot(m1s,scales=list(y=list(rot=0)),main="All Posterior Samples of Model Parameters"))
260
  m1s=mcmc(m1$p.samples,start=istart,thin=ithin)
261
  print(xyplot(m1s,scales=list(y=list(rot=0)),main="Post-burnin, thinned, Posterior Samples of Model Parameters"))
262

  
263

  
264
### Pred vs. obs  plots
265
plot(pd2[,as.character(ivar)],pd2$Q50,pch=16,cex=.8,ylim=quantile(c(pd2$Q2.5,pd2$Q97.5),c(0.01,0.995)),xlim=quantile(c(pd2$Q2.5,pd2$Q97.5),c(0.01,0.995)))
266
segments(pd2[,as.character(ivar)],pd2$Q2.5,pd2[,as.character(ivar)],pd2$Q97.5,col="grey")
267
points(pd2[,as.character(ivar)],pd2$Q50,pch=16,cex=.8);abline(0,1,col="red") #xlim=c(0,10)
268

  
269
### Maps
270
par(mfrow=c(1,3))
271
## observed
272
obs.surf <- mba.surf(cbind(fd[,c("lon","lat")],fd$tmax), no.X=100, no.Y=100, extend=TRUE)$xyz.est
273
image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response",asp=1)
274
points(fd[,c("lon","lat")])
275
plot(roi,add=T)
276
#contour(obs.surf, add=T)
277
plot(roi,add=T,border="darkgreen",lwd=3)
278
## spatial effects
279
  w.pred.surf <-
280
mba.surf(cbind(coordinates(knots), m1p.ws[,3]), no.X=100, no.Y=100, extend=TRUE)$xyz.est
281
image(w.pred.surf, xaxs = "r", yaxs = "r", main="Spatial Effects",asp=1)
282
#points(pd[,c("lon","lat")], pch=1, cex=1)
283
points(knots, pch=3, cex=1)
284
#contour(w.pred.surf, add=T)
285
plot(roi,add=T,border="darkgreen",lwd=3)
286
legend(1.5,2.5, legend=c("Obs.", "Knots", "Pred."),
287
pch=c(1,3,19), bg="white")
288
## fitted
289
y.pred.surf <-
290
mba.surf(cbind(pd2[,c("lon","lat")], pd2$Q50), no.X=100, no.Y=100, extend=TRUE)$xyz.est
291
image(y.pred.surf, xaxs = "r", yaxs = "r", main="Predicted response",asp=1)
292
points(pd[,c("lon","lat")], pch=1, cex=1)
293
#points(knots, pch=3, cex=1)
294
#contour(y.pred.surf, add=T)
295
plot(roi,add=T,border="darkgreen",lwd=3)
296
legend(1.5,2.5, legend=c("Obs.", "Knots", "Pred."),
297
pch=c(1,3,19), bg="white")
298

  
299
dev.off()
300

  
301
  ### Return objects
302
  return(m1.summary)
303
}
304

  
305

  
306

  
307
###########################################################
308
### define models to run
309
ms=expand.grid(variable="tmax",type=c("anom","raw"),date=dates,model=c("intercept","full"))
310
## drop some
311
ms[!ms$type=="raw"&ms$model=="intercept",]
312

  
313

  
314
msl=apply(ms,1,as.list)
315

  
316

  
317
### test it out
318
parms=msl[[9]]
319
interp(parms,niter=5000,istart=2000,ithin=3)
320

  
321
#mresults=lapply(msl[1:3],interp)
322

  
323
### run it!
324
mresults =mclapply(msl,interp,niter=1000,istart=500,ithin=1,mc.cores=1)
325

  
326
### drop any errors
327
mresults=mresults[!sapply(mresults,function(x) grepl("Error",x[[1]]))]
328

  
329
res=do.call(rbind.data.frame,mresults)
330

  
331
save(mresults,res,file="output/mresults.Rdata")
332

  
333
q("no")
334

  
335

  
336

  
337

  
338

  
climate/research/LST_landcover_exploration.R
1
library(sp)                                                                                                                                                                          
2
library(spgrass6)                                                                                                                                                                    
3
library(rgdal)                                                                                                                                                                       
4
library(reshape)                                                                                                                                                                     
5
library(ncdf4)                                                                                                                                                                       
6
library(geosphere)                                                                                                                                                                   
7
library(rgeos)                                                                                                                                                                       
8
library(multicore)                                                                                                                                                                   
9
library(raster)                                                                                                                                                                      
10
library(lattice)                                                                                                                                                                     
11
library(rgl)                                                                                                                                                                         
12
library(hdf5)                                                                                                                                                                        
13
library(rasterVis)                                                                                                                                                                   
14
library(heR.Misc)
15
library(spBayes)
16

  
17

  
18
X11.options(type="X11")
19

  
20
ncores=20  #number of threads to use
21

  
22

  
23
### copy lulc data to litoria
24
setwd("data/lulc")
25
system("scp atlas:/home/parmentier/data_Oregon_stations/W_Layer* .")
26

  
27

  
28
setwd("/home/adamw/acrobates/projects/interp")       
29

  
30
projs=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")
31
months=format(ISOdate(2004,1:12,1),"%b")
32

  
33
### load station data, subset to stations in region, and transform to sinusoidal
34
load("data/ghcn/roi_ghcn.Rdata")
35
load("data/allstations.Rdata")
36

  
37
d2=d[d$variable=="tmax"&d$date>=as.Date("2000-01-01"),]
38
d2=d2[,-grep("variable",colnames(d2)),]
39
st2=st[st$id%in%d$id,]
40
st2=spTransform(st2,projs)
41
d2[,c("lon","lat")]=coordinates(st2)[match(d2$id,st2$id),]
42
d2$elev=st2$elev[match(d2$id,st2$id)]
43
d2$month=format(d2$date,"%m")
44
#d2$value=d2$value/10 #convert to mm
45

  
46

  
47
## load topographical data
48
topo=brick(as.list(list.files("data/topography",pattern="rst$",full=T)))
49
topo=calc(topo,function(x) ifelse(x<0,NA,x))
50
names(topo)=c("aspect","dem","slope")
51
colnames(topo@data@values)=c("aspect","dem","slope")
52
projection(topo)=projs
53
NS=sin(subset(topo,subset="slope")*pi/180)*cos(subset(topo,subset="aspect")*pi/180);names(NS)="northsouth"
54
EW=sin(subset(topo,subset="slope")*pi/180)*sin(subset(topo,subset="aspect")*pi/180);names(EW)="eastwest"
55
slope=sin(subset(topo,subset="slope")*pi/180);names(slope)="slope"
56
topo2=stack(subset(topo,"dem"),EW,NS,slope)
57

  
58
## create binned elevation for stratified sampling
59
demc=quantile(subset(topo,subset="dem"))
60
demb=calc(subset(topo,subset="dem"),function(x) as.numeric(cut(x,breaks=demc)))
61
names(demb)="demb"
62
                                        #topo=brick(list(topo,demb))
63

  
64

  
65
### load the lulc data as a brick
66
lulc=brick(as.list(list.files("data/lulc",pattern="rst$",full=T)))
67
#projection(lulc)=
68
#plot(lulc)
69

  
70
## Enter lulc types (from Nakaegawa 2011)
71
lulct=as.data.frame(matrix(c(
72
  "Forest",1,
73
  "Shrub",2,
74
  "Grass",3,
75
  "Crop",4,
76
  "Mosaic",5,
77
  "Urban",6,
78
  "Barren",7,
79
  "Snow",8,
80
  "Wetland",9,
81
  "Water",10),byrow=T,ncol=2,dimnames=list(1:10,c("class","id"))),stringsAsFactors=F)
82
colnames(lulc@data@values)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
83
layerNames(lulc)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
84
lulc=calc(lulc,function(x) ifelse(is.na(x),0,x))
85
projection(lulc)=projs
86

  
87
### load the LST data
88
lst=brick(as.list(list.files("data/lst",pattern="rescaled.rst$",full=T)[c(4:12,1:3)]))
89
lst=lst-273.15
90
colnames(lst@data@values)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
91
layerNames(lst)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
92
projection(lst)=projs
93

  
94

  
95
######################################
96
## compare LULC with station data
97
stlulc=extract(lulc,st2) #overlay stations and LULC values
98
st2$lulc=do.call(c,lapply(apply(stlulc,1,function(x) which.max(x)),function(x) ifelse(is.null(names(x)),NA,names(x))))
99

  
100

  
101
### generate sample of points to speed processing
102
n=10000/length(unique(demb))
103
n2=30  #number of knots
104
s=sampleStratified(demb,size=n,sp=T)
105
#s=spsample(as(topo,"SpatialGrid"),n=n,type="regular")
106
s2=spsample(as(topo,"SpatialGrid"),n=n2,type="regular")
107

  
108
         
109
s=SpatialPointsDataFrame(s,data=cbind.data.frame(extract(topo,s),extract(topo2,s),extract(lulc,s),extract(lst,s)))
110
### drop areas with no LST data
111
#s=s[!is.na(s$Aug),]
112
s=s[apply(s@data,1,function(x) all(!is.na(x))),]
113

  
114
###  add majority rules lulc for exploration
115
s$lulc=apply(s@data[,colnames(s@data)%in%lulct$class],1,function(x) (colnames(s@data)[colnames(s@data)%in%lulct$class])[which.max(x)])
116

  
117

  
118
### spatial regression to fit lulc coefficients
119

  
120
niter=1250
121
mclapply(months,function(m){
122
  print(paste("###################################   Starting month ",m))
123
  f1=formula(paste(m,"~Shrub+Grass+Crop+Mosaic+Urban+Barren+Snow+Wetland+dem+eastwest+northsouth",sep=""))
124
  m1=spLM(f1,data=s@data,coords=coordinates(s),knots=coordinates(s2),
125
    starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
126
    sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
127
    priors=list("phi.Unif"=c(0.3, 3), "sigma.sq.IG"=c(2, 1), "tau.sq.IG"=c(2, 1)),
128
    cov.model="exponential",
129
    n.samples=niter, verbose=TRUE, n.report=100)
130
  assign(paste("mod_",m,sep=""),m1)
131
  save(list=paste("mod_",m,sep=""),file=paste("output/mod_",m,".Rdata",sep=""))
132
})
133

  
134

  
135
m1.s=mcmc(t(m1$p.samples),start=round(niter/4),thin=1)
136
bwplot(value~X2,melt(as.matrix(m1$p.samples)[,!colnames(m1$p.samples)%in%c("(Intercept)","sigma.sq","tau.sq","phi","northsouth","dem","eastwest")]))
137
densityplot(m1$p.samples)
138

  
139
### load oregon boundary for comparison
140
roi=spTransform(as(readOGR("data/regions/Test_sites/Oregon.shp","Oregon"),"SpatialLines"),projs)
141

  
142

  
143
bgyr=colorRampPalette(c("blue","green","yellow","red"))
144
pdf("output/lst_lulc.pdf",width=11,height=8.5)
145

  
146
### Summary plots of covariates
147
## LULC
148
at=seq(0.1,100,length=100)
149
levelplot(lulc,at=at,col.regions=bgyr(length(at)),
150
          main="Land Cover Classes",sub="Sub-pixel %")+
151
  layer(sp.lines(roi, lwd=1.2, col='black'))
152

  
153
## TOPO
154
at=unique(quantile(as.matrix(subset(topo2,c("eastwest","northsouth","slope"))),seq(0,1,length=100),na.rm=T))
155
levelplot(subset(topo2,c("eastwest","northsouth","slope")),at=at,col.regions=bgyr(length(at)),
156
          main="Topographic Variables",sub="")+
157
  layer(sp.lines(roi, lwd=1.2, col='black'))
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff