kwtrp\lolweny — Sep 28, 2015, 6:36 PM
#Building Spatial Knots from the fields package.
#Borrowing from the philosophy of "PrevMap".
rm(list=ls())
setwd('~/stan_mcmc')
library(geoR)
Warning in fun(libname, pkgname): couldn't connect to display ":0"
--------------------------------------------------------------
Analysis of Geostatistical Data
For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
geoR version 1.7-5.1 (built on 2015-04-15) is now loaded
--------------------------------------------------------------
library(raster)
Loading required package: sp
Warning: no function found corresponding to methods exports from 'raster'
for: 'overlay'
library(sp)
library(maptools)
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()
library(Rcpp)
library(INLA)
Loading required package: Matrix
Loading required package: splines
library(rstan)
Loading required package: inline
Attaching package: 'inline'
The following object is masked from 'package:Rcpp':
registerPlugin
rstan (Version 2.7.0-1, packaged: 2015-07-17 18:12:01 UTC, GitRev: 05c3d0058b6a)
For execution on a local, multicore CPU with excess RAM we recommend calling
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
Attaching package: 'rstan'
The following object is masked from 'package:raster':
extract
library(stringr)
#parallel computing on clusters
library(parallel)
library(Rmpi) #using MPI parallelization scheme
Error in library(Rmpi): there is no package called 'Rmpi'
library("PrevMap") #for low-rank approaximation, for mesh-free methods
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: pdist
library("fields") #for low-rank approaximation
Loading required package: spam
Loading required package: grid
Spam version 1.0-1 (2014-09-09) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: 'spam'
The following objects are masked from 'package:base':
backsolve, forwardsolve
Loading required package: maps
library("maxLik") #for optimizer maxBFGS
library(stats)#for nlminb (Unconstrained optimization )
library(spam)
#package fields
#the 188 exploratory topics
#add.image
#function Krig with A Knot example:
data(ozone2)
y16<- ozone2$y[16,]
good<- !is.na(y16)
y<- y16[good]
x<- ozone2$lon.lat[good,] #147 data points
#
# the knots can be arbitrary but just for fun find them with a space
# filling design. Here we select 50 from the full set of 147 points
#
#knots
#one way may be space-filling the space spanned by the observations
#or putting knots at all inique locations
knots<- cover.design( x, 50, num.nn= 75)
xknots<- cover.design( x, 50, num.nn= 75)$design # select 50 knot points/50 basis functions
out<- Krig( x, y, knots=xknots, cov.function="Exp.cov", theta=300)
summary( out)
CALL:
Krig(x = x, Y = y, cov.function = "Exp.cov", knots = xknots,
theta = 300)
Number of Observations: 147
Number of unique points: 147
Number of parameters in the null space 3
Parameters for fixed spatial drift 3
Effective degrees of freedom: 34.8
Residual degrees of freedom: 112.2
MLE sigma NA
GCV sigma 12.49
MLE rho NA
Scale passed for covariance (rho) <NA>
Scale passed for nugget (sigma^2) <NA>
Smoothing parameter lambda 0.002346
Residual Summary:
min 1st Q median 3rd Q max
-39.550 -6.207 0.449 5.965 31.800
Covariance Model: Exp.cov
Names of non-default covariance arguments:
theta
Knot model: 50 knots supplied to define basis
functions
DETAILS ON SMOOTHING PARAMETER:
Method used: REML Cost: 1
lambda trA GCV GCV.one GCV.model shat
2.346e-03 3.485e+01 1.776e+03 2.045e+02 NA 1.249e+01
Summary of all estimates found for lambda
lambda trA GCV shat -lnLike Prof converge
GCV 0.019444 16.47 897.9 14.92 329.1 8
GCV.model NA NA NA NA NA NA
GCV.one 0.002826 33.22 204.2 12.57 NA 6
RMSE NA NA NA NA NA NA
pure error NA NA NA NA NA NA
REML 0.002346 34.85 1775.5 12.49 325.2 3
cbind(out$fitted.values,y)
y
[1,] 76.844096 75.00000
[2,] 85.165805 84.25000
[3,] 95.170347 90.87500
[4,] 105.394517 127.42857
[5,] 93.500201 104.50000
[6,] 94.362577 86.25000
[7,] 90.623601 93.25000
[8,] 105.091485 95.50000
[9,] 88.971206 89.87500
[10,] 98.913241 100.87500
[11,] 106.794081 87.25000
[12,] 110.000748 112.37500
[13,] 89.625164 91.37500
[14,] 82.434976 88.25000
[15,] 98.896146 89.12500
[16,] 116.440957 115.12500
[17,] 132.093586 140.87500
[18,] 120.292316 115.75000
[19,] 106.519424 98.00000
[20,] 87.614209 87.37500
[21,] 89.860748 89.37500
[22,] 78.670875 82.87500
[23,] 77.028061 76.75000
[24,] 77.648436 78.75000
[25,] 78.426023 78.87500
[26,] 70.891556 59.00000
[27,] 92.211663 87.75000
[28,] 91.742953 100.75000
[29,] 82.355168 76.00000
[30,] 73.773764 58.00000
[31,] 90.608607 89.50000
[32,] 87.506315 91.00000
[33,] 83.322911 83.50000
[34,] 83.371289 83.50000
[35,] 71.727404 60.50000
[36,] 70.772433 73.62500
[37,] 73.855064 73.75000
[38,] 63.656301 48.00000
[39,] 87.643184 78.85714
[40,] 88.252365 86.50000
[41,] 88.887020 88.33333
[42,] 66.454180 76.00000
[43,] 71.626739 48.00000
[44,] 88.192792 67.87500
[45,] 89.155272 92.00000
[46,] 90.219477 97.00000
[47,] 89.418658 88.71429
[48,] 89.904694 96.87500
[49,] 89.169719 95.00000
[50,] 87.379087 95.37500
[51,] 86.454031 76.87500
[52,] 85.006698 71.25000
[53,] 85.373573 98.25000
[54,] 78.838820 72.62500
[55,] 57.092765 65.25000
[56,] 59.144739 67.00000
[57,] 87.034846 96.28571
[58,] 58.721505 69.85714
[59,] 65.117582 79.66667
[60,] 63.985875 55.75000
[61,] 65.241918 66.50000
[62,] 8.537520 0.00000
[63,] 9.455186 5.25000
[64,] 82.671200 87.62500
[65,] 86.300677 96.87500
[66,] 54.815747 55.75000
[67,] 94.560384 89.62500
[68,] 59.828547 61.25000
[69,] 43.751042 40.85714
[70,] 68.907398 71.87500
[71,] 67.716439 59.37500
[72,] 60.369967 72.75000
[73,] 55.690778 22.75000
[74,] 56.236674 53.87500
[75,] 60.390580 60.62500
[76,] 94.294427 100.50000
[77,] 53.134607 58.12500
[78,] 52.227161 57.50000
[79,] 40.900132 35.25000
[80,] 65.344431 68.25000
[81,] 38.087941 30.37500
[82,] 78.999873 85.33333
[83,] 72.401987 73.00000
[84,] 70.154903 68.40000
[85,] 80.576948 70.40000
[86,] 97.822410 100.75000
[87,] 96.175338 100.50000
[88,] 61.473440 66.75000
[89,] 107.055001 96.37500
[90,] 104.266989 103.16667
[91,] 59.955613 63.75000
[92,] 77.084532 95.75000
[93,] 75.094616 78.37500
[94,] 63.786458 52.87500
[95,] 61.421277 51.62500
[96,] 61.467637 62.75000
[97,] 66.753557 54.75000
[98,] 55.191873 58.28571
[99,] 77.899946 84.00000
[100,] 72.337551 79.85714
[101,] 73.145821 82.00000
[102,] 73.763529 82.50000
[103,] 75.791813 91.25000
[104,] 74.914202 95.87500
[105,] 72.668557 60.62500
[106,] 73.540664 65.25000
[107,] 73.598304 60.00000
[108,] 74.516313 68.37500
[109,] 85.338489 93.00000
[110,] 84.022707 81.50000
[111,] 90.230851 89.14286
[112,] 88.200798 87.87500
[113,] 97.238719 106.12500
[114,] 96.752779 84.25000
[115,] 95.153315 108.37500
[116,] 92.448317 86.57143
[117,] 92.967217 96.75000
[118,] 96.231206 106.50000
[119,] 66.215043 62.57143
[120,] 66.620877 64.62500
[121,] 82.174181 73.50000
[122,] 85.902016 86.28571
[123,] 89.138639 88.87500
[124,] 79.002179 80.87500
[125,] 81.930152 84.37500
[126,] 82.394014 82.87500
[127,] 90.752784 96.00000
[128,] 93.477811 80.37500
[129,] 87.578698 89.62500
[130,] 130.772769 162.57143
[131,] 118.824640 112.62500
[132,] 63.631790 66.62500
[133,] 128.116358 137.75000
[134,] 128.539675 157.62500
[135,] 126.882165 87.33333
[136,] 127.418050 107.75000
[137,] 126.964546 157.37500
[138,] 128.957960 137.37500
[139,] 90.303346 72.50000
[140,] 123.374797 148.50000
[141,] 130.201671 144.00000
[142,] 84.240121 82.87500
[143,] 86.494310 84.87500
[144,] 124.439546 144.62500
[145,] 108.866097 94.75000
[146,] 116.055032 88.00000
[147,] 81.116599 81.62500
#extracting predictions from out
out$m
[1] 2
out$chol.args
$pivot
[1] FALSE
out$beta
NULL
out$transform
$dim
[1] 147 2
$`scaled:center`
[1] 0 0
$`scaled:scale`
[1] 1 1
$x.center
[1] 0 0
$x.scale
[1] 1 1
$x.scale.type
[1] "user"
out$gcv.grid
lambda trA GCV GCV.one GCV.model shat
1 3.052855e-05 50.350007 66025.8988 239.4041 NA 12.54608
2 4.058313e-05 50.112057 55867.4678 238.9040 NA 12.54839
3 5.253526e-05 49.874127 47909.4318 238.3760 NA 12.54989
4 6.626474e-05 49.636172 41546.7419 237.7941 NA 12.54991
5 8.155845e-05 49.398262 36375.3429 237.1499 NA 12.54820
6 9.819145e-05 49.160311 32111.7430 236.4470 NA 12.54486
7 1.159531e-04 48.922350 28555.4174 235.6953 NA 12.54012
8 1.346743e-04 48.684367 25558.4633 234.9059 NA 12.53429
9 1.542117e-04 48.446516 23011.5726 234.0900 NA 12.52763
10 1.744888e-04 48.208548 20827.3383 233.2556 NA 12.52037
11 1.954215e-04 47.970622 18941.7134 232.4103 NA 12.51270
12 2.169657e-04 47.732659 17302.4921 231.5597 NA 12.50478
13 2.390784e-04 47.494730 15869.2391 230.7083 NA 12.49672
14 2.617400e-04 47.256770 14608.7128 229.8594 NA 12.48861
15 2.849269e-04 47.018838 13494.6770 229.0161 NA 12.48054
16 3.086297e-04 46.780911 12505.3115 228.1802 NA 12.47256
17 3.328423e-04 46.542975 11622.7452 227.3535 NA 12.46471
18 3.575598e-04 46.305044 10832.2658 226.5373 NA 12.45705
19 3.827837e-04 46.067094 10121.4961 225.7327 NA 12.44959
20 4.085123e-04 45.829152 9480.1965 224.9404 NA 12.44236
21 4.347503e-04 45.591202 8899.6098 224.1612 NA 12.43539
22 4.615033e-04 45.353234 8372.3265 223.3956 NA 12.42869
23 4.887663e-04 45.115341 7892.2315 222.6442 NA 12.42228
24 5.165630e-04 44.877372 7453.6190 221.9070 NA 12.41618
25 5.448893e-04 44.639431 7052.0401 221.1845 NA 12.41038
26 5.737530e-04 44.401517 6683.4644 220.4769 NA 12.40490
27 6.031617e-04 44.163643 6344.4071 219.7845 NA 12.39976
28 6.331467e-04 43.925635 6031.5975 219.1068 NA 12.39494
29 6.636937e-04 43.687689 5742.6456 218.4445 NA 12.39047
30 6.948308e-04 43.449673 5475.0561 217.7973 NA 12.38635
31 7.265405e-04 43.211803 5227.0082 217.1658 NA 12.38258
32 7.588602e-04 42.973886 4996.4740 216.5496 NA 12.37916
33 7.918007e-04 42.735937 4781.8680 215.9486 NA 12.37610
34 8.253702e-04 42.497988 4581.7960 215.3630 NA 12.37341
35 8.595819e-04 42.260042 4394.9869 214.7928 NA 12.37107
36 8.944540e-04 42.022069 4220.2832 214.2379 NA 12.36911
37 9.299787e-04 41.784217 4056.7693 213.6986 NA 12.36752
38 9.662052e-04 41.546256 3903.3741 213.1743 NA 12.36630
39 1.003136e-03 41.308271 3759.3396 212.6653 NA 12.36546
40 1.040775e-03 41.070343 3623.9725 212.1715 NA 12.36499
41 1.079145e-03 40.832429 3496.5759 211.6930 NA 12.36490
42 1.118272e-03 40.594474 3376.5174 211.2296 NA 12.36519
43 1.158164e-03 40.356530 3263.2752 210.7812 NA 12.36586
44 1.198841e-03 40.118588 3156.3465 210.3479 NA 12.36692
45 1.240319e-03 39.880661 3055.2826 209.9297 NA 12.36836
46 1.282633e-03 39.642667 2959.6366 209.5263 NA 12.37019
47 1.325774e-03 39.404758 2869.0912 209.1380 NA 12.37241
48 1.369780e-03 39.166840 2783.2631 208.7646 NA 12.37502
49 1.414666e-03 38.928939 2701.8450 208.4060 NA 12.37802
50 1.460466e-03 38.690999 2624.5277 208.0622 NA 12.38141
51 1.507202e-03 38.453019 2551.0466 207.7330 NA 12.38520
52 1.554883e-03 38.215070 2481.1776 207.4187 NA 12.38938
53 1.603538e-03 37.977131 2414.6879 207.1191 NA 12.39396
54 1.653212e-03 37.739098 2351.3411 206.8341 NA 12.39895
55 1.703849e-03 37.501357 2291.0440 206.5641 NA 12.40432
56 1.755598e-03 37.263327 2233.4706 206.3083 NA 12.41011
57 1.808410e-03 37.025373 2178.5480 206.0673 NA 12.41630
58 1.862275e-03 36.787652 2126.1544 205.8410 NA 12.42288
59 1.917369e-03 36.549519 2076.0089 205.6290 NA 12.42989
60 1.973587e-03 36.311563 2028.1068 205.4317 NA 12.43730
61 2.031024e-03 36.073513 1982.2696 205.2489 NA 12.44513
62 2.089637e-03 35.835672 1938.4403 205.0807 NA 12.45336
63 2.149526e-03 35.597766 1896.4601 204.9271 NA 12.46201
64 2.210725e-03 35.359793 1856.2312 204.7880 NA 12.47108
65 2.273258e-03 35.121806 1817.6710 204.6635 NA 12.48056
66 2.337130e-03 34.883915 1780.7096 204.5535 NA 12.49047
67 2.402418e-03 34.645976 1745.2435 204.4581 NA 12.50080
68 2.469144e-03 34.408048 1711.2063 204.3774 NA 12.51156
69 2.537359e-03 34.170097 1678.5221 204.3112 NA 12.52275
70 2.607094e-03 33.932155 1647.1292 204.2597 NA 12.53436
71 2.678385e-03 33.694248 1616.9681 204.2229 NA 12.54641
72 2.751341e-03 33.456164 1587.9546 204.2008 NA 12.55891
73 2.825868e-03 33.218355 1560.0880 204.1935 NA 12.57183
74 2.902141e-03 32.980415 1533.2680 204.2010 NA 12.58519
75 2.980212e-03 32.742344 1507.4476 204.2233 NA 12.59902
76 3.060013e-03 32.504498 1482.6193 204.2605 NA 12.61327
77 3.141695e-03 32.266586 1458.7084 204.3127 NA 12.62798
78 3.225304e-03 32.028636 1435.6776 204.3800 NA 12.64315
79 3.310875e-03 31.790710 1413.4947 204.4623 NA 12.65878
80 3.398505e-03 31.552704 1392.1142 204.5599 NA 12.67487
81 3.488170e-03 31.314853 1371.5224 204.6727 NA 12.69142
82 3.580047e-03 31.076854 1351.6610 204.8008 NA 12.70844
83 3.674100e-03 30.838971 1332.5214 204.9444 NA 12.72593
84 3.770482e-03 30.600997 1314.0580 205.1036 NA 12.74391
85 3.869187e-03 30.363117 1296.2579 205.2783 NA 12.76235
86 3.970385e-03 30.125104 1279.0783 205.4689 NA 12.78130
87 4.074061e-03 29.887180 1262.5110 205.6753 NA 12.80073
88 4.180333e-03 29.649252 1246.5258 205.8977 NA 12.82065
89 4.289284e-03 29.411324 1231.1011 206.1362 NA 12.84107
90 4.401022e-03 29.173351 1216.2133 206.3911 NA 12.86200
91 4.515596e-03 28.935422 1201.8481 206.6623 NA 12.88344
92 4.633147e-03 28.697442 1187.9809 206.9501 NA 12.90539
93 4.753698e-03 28.459560 1174.6024 207.2545 NA 12.92786
94 4.877431e-03 28.221619 1161.6868 207.5758 NA 12.95085
95 5.001045e-03 27.989928 1149.5432 207.9051 NA 12.97376
96 5.134793e-03 27.745725 1137.1902 208.2699 NA 12.99845
97 5.268637e-03 27.507794 1125.5804 208.6429 NA 13.02306
98 5.406084e-03 27.269871 1114.3774 209.0336 NA 13.04821
99 5.547285e-03 27.031915 1103.5661 209.4421 NA 13.07393
100 5.692341e-03 26.793980 1093.1366 209.8686 NA 13.10021
101 5.841417e-03 26.556021 1083.0750 210.3134 NA 13.12706
102 5.994619e-03 26.318099 1073.3726 210.7766 NA 13.15448
103 6.152075e-03 26.080246 1064.0195 211.2584 NA 13.18247
104 6.314142e-03 25.842166 1054.9944 211.7596 NA 13.21109
105 6.480710e-03 25.604267 1046.3029 212.2797 NA 13.24028
106 6.652104e-03 25.366326 1037.9275 212.8193 NA 13.27009
107 6.828454e-03 25.128414 1029.8617 213.3787 NA 13.30050
108 7.010017e-03 24.890442 1022.0943 213.9582 NA 13.33155
109 7.196912e-03 24.652516 1014.6207 214.5580 NA 13.36322
110 7.389390e-03 24.414579 1007.4316 215.1785 NA 13.39554
111 7.587669e-03 24.176636 1000.5198 215.8200 NA 13.42850
112 7.791963e-03 23.938703 993.8789 216.4827 NA 13.46213
113 8.002540e-03 23.700751 987.5017 217.1672 NA 13.49642
114 8.219667e-03 23.462767 981.3817 217.8739 NA 13.53140
115 8.443488e-03 23.224883 975.5162 218.6026 NA 13.56706
116 8.674482e-03 22.986889 969.8946 219.3545 NA 13.60343
117 8.912726e-03 22.749008 964.5169 220.1292 NA 13.64049
118 9.158701e-03 22.511072 959.3741 220.9277 NA 13.67829
119 9.412695e-03 22.273118 954.4625 221.7503 NA 13.71682
120 9.675032e-03 22.035170 949.7780 222.5973 NA 13.75610
121 9.946022e-03 21.797275 945.3174 223.4692 NA 13.79612
122 1.022620e-02 21.559298 941.0742 224.3667 NA 13.83693
123 1.051584e-02 21.321362 937.0469 225.2901 NA 13.87852
124 1.081541e-02 21.083427 933.2313 226.2399 NA 13.92090
125 1.112539e-02 20.845475 929.6238 227.2168 NA 13.96410
126 1.144621e-02 20.607550 926.2222 228.2212 NA 14.00812
127 1.177848e-02 20.369576 923.0225 229.2539 NA 14.05299
128 1.212256e-02 20.131670 920.0237 230.3152 NA 14.09870
129 1.247920e-02 19.893729 917.2221 231.4060 NA 14.14530
130 1.284895e-02 19.655773 914.6159 232.5269 NA 14.19278
131 1.323243e-02 19.417835 912.2033 233.6786 NA 14.24117
132 1.363039e-02 19.179864 909.9823 234.8619 NA 14.29049
133 1.404341e-02 18.941959 907.9521 236.0772 NA 14.34075
134 1.447250e-02 18.703980 906.1103 237.3259 NA 14.39198
135 1.491819e-02 18.466094 904.4571 238.6080 NA 14.44417
136 1.538134e-02 18.228298 902.9912 239.9243 NA 14.49735
137 1.586371e-02 17.990196 901.7099 241.2780 NA 14.55163
138 1.636560e-02 17.752144 900.6148 242.6682 NA 14.60695
139 1.688770e-02 17.514318 899.7064 244.0946 NA 14.66328
140 1.743180e-02 17.276411 898.9832 245.5601 NA 14.72074
141 1.799856e-02 17.038676 898.4462 247.0643 NA 14.77929
142 1.859078e-02 16.800508 898.0947 248.6123 NA 14.83909
143 1.920832e-02 16.562555 897.9308 250.2011 NA 14.90003
144 1.985315e-02 16.324635 897.9550 251.8331 NA 14.96217
145 2.052710e-02 16.086674 898.1686 253.5103 NA 15.02557
146 2.123179e-02 15.848738 898.5731 255.2336 NA 15.09025
147 2.196912e-02 15.610828 899.1703 257.0045 NA 15.15624
148 2.274159e-02 15.372812 899.9628 258.8256 NA 15.22361
149 2.355085e-02 15.134879 900.9524 260.6972 NA 15.29236
150 2.439939e-02 14.896997 902.1419 262.6213 NA 15.36253
151 2.529032e-02 14.659042 903.5352 264.6009 NA 15.43420
152 2.622626e-02 14.421084 905.1355 266.6374 NA 15.50740
153 2.721014e-02 14.183171 906.9462 268.7324 NA 15.58217
154 2.824560e-02 13.945242 908.9722 270.8890 NA 15.65857
155 2.933640e-02 13.707293 911.2184 273.1094 NA 15.73667
156 3.048660e-02 13.469326 913.6899 275.3963 NA 15.81652
157 3.170045e-02 13.231389 916.3920 277.7521 NA 15.89817
158 3.298276e-02 12.993490 919.3308 280.1794 NA 15.98168
159 3.433948e-02 12.755531 922.5143 282.6825 NA 16.06715
160 3.577626e-02 12.517567 925.9493 285.2640 NA 16.15465
161 3.729934e-02 12.279652 929.6428 287.9270 NA 16.24423
162 3.891649e-02 12.041707 933.6049 290.6761 NA 16.33600
163 4.063542e-02 11.803789 937.8438 293.5147 NA 16.43004
164 4.246562e-02 11.565832 942.3710 296.4478 NA 16.52645
165 4.441673e-02 11.327894 947.1963 299.4794 NA 16.62532
166 4.650090e-02 11.089874 952.3336 302.6158 NA 16.72681
167 4.872874e-02 10.851999 957.7907 305.8591 NA 16.83091
168 5.111631e-02 10.614070 963.5855 309.2177 NA 16.93785
169 5.367927e-02 10.376139 969.7318 312.6968 NA 17.04772
170 5.643596e-02 10.138208 976.2454 316.3030 NA 17.16067
171 5.940744e-02 9.900260 983.1438 320.0435 NA 17.27683
172 6.261761e-02 9.662301 990.4455 323.9256 NA 17.39638
173 6.609324e-02 9.424381 998.1686 327.9565 NA 17.51944
174 6.986701e-02 9.186435 1006.3368 332.1458 NA 17.64622
175 7.397547e-02 8.948489 1014.9723 336.5026 NA 17.77690
176 7.888733e-02 8.689042 1024.9499 341.4553 NA 17.92407
177 8.338280e-02 8.472297 1033.7593 345.7646 NA 18.05094
178 8.877881e-02 8.234669 1043.9413 350.6799 NA 18.19438
179 9.474095e-02 7.996740 1054.7156 355.8134 NA 18.34277
180 1.013484e-01 7.758811 1066.1040 361.1727 NA 18.49621
181 1.087060e-01 7.520851 1078.1451 366.7733 NA 18.65498
182 1.169390e-01 7.282920 1090.8765 372.6298 NA 18.81936
183 1.262059e-01 7.044968 1104.3447 378.7608 NA 18.98970
184 1.367011e-01 6.807045 1118.5950 385.1842 NA 19.16632
185 1.486635e-01 6.569325 1133.6680 391.9153 NA 19.34944
186 1.624551e-01 6.331047 1149.6716 398.9993 NA 19.54009
187 1.784410e-01 6.093169 1166.6054 406.4330 NA 19.73794
188 1.972040e-01 5.855253 1184.5678 414.2567 NA 19.94382
189 2.195031e-01 5.617332 1203.6327 422.4993 NA 20.15823
190 2.464049e-01 5.379402 1223.8847 431.1944 NA 20.38173
191 2.794496e-01 5.141467 1245.4168 440.3785 NA 20.61494
192 3.209550e-01 4.903487 1268.3355 450.0940 NA 20.85858
193 3.745279e-01 4.665587 1292.7425 460.3805 NA 21.11323
194 4.462494e-01 4.427627 1318.7778 471.2939 NA 21.37986
195 5.469983e-01 4.189677 1346.5742 482.8866 NA 21.65926
196 6.980708e-01 3.952348 1376.2072 495.1871 NA 21.95161
197 9.517268e-01 3.713826 1408.0823 508.3606 NA 22.26022
198 1.459127e+00 3.475861 1442.1640 522.3892 NA 22.58400
199 2.982912e+00 3.237947 1478.7262 537.3835 NA 22.92480
200 7.299828e+02 3.000994 1517.8434 553.3721 NA 23.28250
-lnLike Prof
1 358.5222
2 354.4339
3 350.8122
4 347.6521
5 344.9305
6 342.6034
7 340.6183
8 338.9225
9 337.4687
10 336.2144
11 335.1266
12 334.1769
13 333.3433
14 332.6071
15 331.9536
16 331.3707
17 330.8483
18 330.3782
19 329.9536
20 329.5686
21 329.2184
22 328.8990
23 328.6070
24 328.3391
25 328.0930
26 327.8664
27 327.6574
28 327.4641
29 327.2854
30 327.1196
31 326.9660
32 326.8232
33 326.6905
34 326.5670
35 326.4521
36 326.3450
37 326.2454
38 326.1524
39 326.0658
40 325.9851
41 325.9099
42 325.8399
43 325.7747
44 325.7141
45 325.6578
46 325.6055
47 325.5571
48 325.5123
49 325.4709
50 325.4328
51 325.3978
52 325.3658
53 325.3366
54 325.3101
55 325.2861
56 325.2647
57 325.2456
58 325.2289
59 325.2143
60 325.2018
61 325.1914
62 325.1830
63 325.1765
64 325.1718
65 325.1690
66 325.1679
67 325.1684
68 325.1707
69 325.1746
70 325.1800
71 325.1870
72 325.1954
73 325.2054
74 325.2168
75 325.2296
76 325.2438
77 325.2594
78 325.2763
79 325.2945
80 325.3141
81 325.3349
82 325.3571
83 325.3804
84 325.4051
85 325.4309
86 325.4580
87 325.4863
88 325.5158
89 325.5465
90 325.5784
91 325.6115
92 325.6457
93 325.6812
94 325.7177
95 325.7545
96 325.7944
97 325.8345
98 325.8757
99 325.9181
100 325.9617
101 326.0064
102 326.0523
103 326.0993
104 326.1475
105 326.1969
106 326.2474
107 326.2992
108 326.3521
109 326.4062
110 326.4615
111 326.5179
112 326.5756
113 326.6345
114 326.6947
115 326.7560
116 326.8186
117 326.8825
118 326.9476
119 327.0139
120 327.0816
121 327.1505
122 327.2208
123 327.2924
124 327.3653
125 327.4395
126 327.5152
127 327.5922
128 327.6706
129 327.7504
130 327.8316
131 327.9143
132 327.9984
133 328.0841
134 328.1712
135 328.2599
136 328.3500
137 328.4419
138 328.5354
139 328.6303
140 328.7270
141 328.8252
142 328.9254
143 329.0271
144 329.1306
145 329.2360
146 329.3431
147 329.4520
148 329.5629
149 329.6757
150 329.7903
151 329.9070
152 330.0257
153 330.1464
154 330.2691
155 330.3940
156 330.5211
157 330.6503
158 330.7817
159 330.9153
160 331.0513
161 331.1895
162 331.3301
163 331.4731
164 331.6184
165 331.7662
166 331.9164
167 332.0690
168 332.2241
169 332.3817
170 332.5418
171 332.7043
172 332.8693
173 333.0367
174 333.2065
175 333.3787
176 333.5689
177 333.7299
178 333.9083
179 334.0889
180 334.2712
181 334.4551
182 334.6402
183 334.8262
184 335.0126
185 335.1988
186 335.3848
187 335.5688
188 335.7503
189 335.9276
190 336.0991
191 336.2623
192 336.4137
193 336.5485
194 336.6596
195 336.7358
196 336.7583
197 336.6917
198 336.4562
199 335.7972
200 327.9673
out$fitted.values.null
[,1]
[1,] 21.66306
[2,] 35.48052
[3,] 39.77872
[4,] 39.42413
[5,] 39.64340
[6,] 38.98918
[7,] 39.61011
[8,] 38.90840
[9,] 37.79418
[10,] 38.96922
[11,] 38.49786
[12,] 39.47230
[13,] 37.55728
[14,] 33.23088
[15,] 36.89321
[16,] 38.75435
[17,] 39.00934
[18,] 38.34460
[19,] 37.17277
[20,] 32.27215
[21,] 28.07706
[22,] 26.20290
[23,] 26.88020
[24,] 26.58720
[25,] 26.35959
[26,] 25.68671
[27,] 30.01830
[28,] 30.16109
[29,] 26.80297
[30,] 25.92064
[31,] 29.32675
[32,] 37.47723
[33,] 33.63071
[34,] 33.84183
[35,] 50.21226
[36,] 49.74708
[37,] 48.84777
[38,] 45.05577
[39,] 45.01640
[40,] 45.60457
[41,] 44.57574
[42,] 37.31892
[43,] 37.53927
[44,] 40.55650
[45,] 43.60986
[46,] 44.77038
[47,] 43.96559
[48,] 44.52711
[49,] 42.03316
[50,] 41.64946
[51,] 41.15587
[52,] 45.95079
[53,] 45.36226
[54,] 41.70148
[55,] 36.57619
[56,] 36.79966
[57,] 38.51878
[58,] 37.60818
[59,] 22.27402
[60,] 22.26085
[61,] 22.84036
[62,] 14.25097
[63,] 13.72279
[64,] 27.27297
[65,] 49.01267
[66,] 44.58727
[67,] 50.78754
[68,] 38.50894
[69,] 42.34795
[70,] 50.01677
[71,] 49.82415
[72,] 39.40006
[73,] 36.41334
[74,] 45.24481
[75,] 44.77351
[76,] 50.55111
[77,] 32.31609
[78,] 31.46531
[79,] 38.61042
[80,] 46.06064
[81,] 33.63221
[82,] 54.20664
[83,] 57.44176
[84,] 58.52435
[85,] 53.54578
[86,] 48.86201
[87,] 50.02388
[88,] 59.94334
[89,] 46.74318
[90,] 46.84003
[91,] 59.14016
[92,] 56.72776
[93,] 56.32495
[94,] 58.82985
[95,] 59.42315
[96,] 59.89683
[97,] 58.13119
[98,] 11.33029
[99,] 25.85313
[100,] 25.05416
[101,] 24.47636
[102,] 25.14510
[103,] 25.50617
[104,] 25.07443
[105,] 25.41585
[106,] 25.39449
[107,] 25.76276
[108,] 25.60515
[109,] 54.42536
[110,] 53.47882
[111,] 52.05714
[112,] 54.02235
[113,] 57.70222
[114,] 57.88307
[115,] 58.12877
[116,] 51.40790
[117,] 49.91830
[118,] 50.51798
[119,] 56.96995
[120,] 57.26123
[121,] 52.58760
[122,] 50.34226
[123,] 52.23213
[124,] 33.23914
[125,] 33.08847
[126,] 33.31598
[127,] 36.56847
[128,] 37.60699
[129,] 35.36259
[130,] 39.23908
[131,] 42.13572
[132,] 25.63403
[133,] 41.15927
[134,] 39.42101
[135,] 38.90745
[136,] 38.99181
[137,] 39.43231
[138,] 39.29486
[139,] 38.28506
[140,] 39.37631
[141,] 39.48385
[142,] 33.87776
[143,] 34.55459
[144,] 40.52999
[145,] 37.91020
[146,] 37.86105
[147,] 37.47223
out$args
$theta
[1] 300
out$lambda
[1] 0.002346449
out$fixed.model
[1] FALSE
out$fitted.values
[,1]
[1,] 76.844096
[2,] 85.165805
[3,] 95.170347
[4,] 105.394517
[5,] 93.500201
[6,] 94.362577
[7,] 90.623601
[8,] 105.091485
[9,] 88.971206
[10,] 98.913241
[11,] 106.794081
[12,] 110.000748
[13,] 89.625164
[14,] 82.434976
[15,] 98.896146
[16,] 116.440957
[17,] 132.093586
[18,] 120.292316
[19,] 106.519424
[20,] 87.614209
[21,] 89.860748
[22,] 78.670875
[23,] 77.028061
[24,] 77.648436
[25,] 78.426023
[26,] 70.891556
[27,] 92.211663
[28,] 91.742953
[29,] 82.355168
[30,] 73.773764
[31,] 90.608607
[32,] 87.506315
[33,] 83.322911
[34,] 83.371289
[35,] 71.727404
[36,] 70.772433
[37,] 73.855064
[38,] 63.656301
[39,] 87.643184
[40,] 88.252365
[41,] 88.887020
[42,] 66.454180
[43,] 71.626739
[44,] 88.192792
[45,] 89.155272
[46,] 90.219477
[47,] 89.418658
[48,] 89.904694
[49,] 89.169719
[50,] 87.379087
[51,] 86.454031
[52,] 85.006698
[53,] 85.373573
[54,] 78.838820
[55,] 57.092765
[56,] 59.144739
[57,] 87.034846
[58,] 58.721505
[59,] 65.117582
[60,] 63.985875
[61,] 65.241918
[62,] 8.537520
[63,] 9.455186
[64,] 82.671200
[65,] 86.300677
[66,] 54.815747
[67,] 94.560384
[68,] 59.828547
[69,] 43.751042
[70,] 68.907398
[71,] 67.716439
[72,] 60.369967
[73,] 55.690778
[74,] 56.236674
[75,] 60.390580
[76,] 94.294427
[77,] 53.134607
[78,] 52.227161
[79,] 40.900132
[80,] 65.344431
[81,] 38.087941
[82,] 78.999873
[83,] 72.401987
[84,] 70.154903
[85,] 80.576948
[86,] 97.822410
[87,] 96.175338
[88,] 61.473440
[89,] 107.055001
[90,] 104.266989
[91,] 59.955613
[92,] 77.084532
[93,] 75.094616
[94,] 63.786458
[95,] 61.421277
[96,] 61.467637
[97,] 66.753557
[98,] 55.191873
[99,] 77.899946
[100,] 72.337551
[101,] 73.145821
[102,] 73.763529
[103,] 75.791813
[104,] 74.914202
[105,] 72.668557
[106,] 73.540664
[107,] 73.598304
[108,] 74.516313
[109,] 85.338489
[110,] 84.022707
[111,] 90.230851
[112,] 88.200798
[113,] 97.238719
[114,] 96.752779
[115,] 95.153315
[116,] 92.448317
[117,] 92.967217
[118,] 96.231206
[119,] 66.215043
[120,] 66.620877
[121,] 82.174181
[122,] 85.902016
[123,] 89.138639
[124,] 79.002179
[125,] 81.930152
[126,] 82.394014
[127,] 90.752784
[128,] 93.477811
[129,] 87.578698
[130,] 130.772769
[131,] 118.824640
[132,] 63.631790
[133,] 128.116358
[134,] 128.539675
[135,] 126.882165
[136,] 127.418050
[137,] 126.964546
[138,] 128.957960
[139,] 90.303346
[140,] 123.374797
[141,] 130.201671
[142,] 84.240121
[143,] 86.494310
[144,] 124.439546
[145,] 108.866097
[146,] 116.055032
[147,] 81.116599
# note that that trA found by GCV is around 17 so 50>17 knots may be a
# reasonable approximation to the full estimator.
#
## Not run:
# the plot
surface( out, type="C")
US( add=TRUE)
points( x, col=2)
points( xknots, cex=2, pch="O")
## End(Not run)
## Not run:
## A quick way to deal with too much data if you intend to smooth it anyway
## Discretize the locations to a grid, then apply Krig
## to the discretized locations:
##
RM.approx<- as.image(RMprecip$y, x=RMprecip$x, nx=20, ny=20)
# take a look:
image.plot( RM.approx)
# discretized data (observations averaged if in the same grid box)
# 336 locations -- down form the full 806
# convert the image format to locations, obs and weight vectors
yd<- RM.approx$z[RM.approx$ind]
weights<- RM.approx$weights[RM.approx$ind] # takes into account averaging
xd<- RM.approx$xd
obj<- Krig( xd, yd, weights=weights, theta=4)
# compare to the full fit:
# Krig( RMprecip$x, RMprecip$y, theta=4)
## End(Not run)
## Not run:
# A correlation model example
# fit krig surface using a mean and sd function to standardize
# first get stats from 1987 summer Midwest O3 data set
data(ozone2)
stats.o3<- stats( ozone2$y)
mean.o3<- Tps( ozone2$lon.lat, c( stats.o3[2,]))
sd.o3<- Tps( ozone2$lon.lat, c( stats.o3[3,]))
#
# Now use these to fit particular day ( day 16)
# and use great circle distance
fit<- Krig( ozone2$lon.lat, ozone2$y[16,],
theta=350, mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern", Distance="rdist.earth",
smoothness=1.0,
na.rm=TRUE) #
# the finale
surface( fit, type="I")
US( add=TRUE)
points( fit$x)
title("Estimated ozone surface")
## End(Not run)
## Not run:
#
#
# explore some different values for the range and lambda using REML
theta <- seq( 100,500,,40)
PLL<- matrix( NA, 40,80)
# the loop
for( k in 1:40){
# call to Krig with different ranges
# also turn off warnings for GCV search
# to avoid lots of messages. (not recommended in general!)
PLL[k,]<- Krig( ozone2$lon.lat,ozone2$y[16,],
cov.function="stationary.cov",
theta=theta[k], mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern",smoothness=.5,
Distance="rdist.earth", nstep.cv=80,
give.warnings=FALSE, na.rm=TRUE)$gcv.grid[,7]
#
# gcv.grid is the grid search output from
# the optimization for estimating different estimates for lambda including
# REML
# default grid is equally spaced in eff.df scale ( and should the same across theta)
# here
}
# get lambda grid from looping
k<- 1
lam<- Krig( ozone2$lon.lat,ozone2$y[16,],
cov.function="stationary.cov",
theta=theta[k], mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern",smoothness=.5,
Distance="rdist.earth", nstep.cv=80,
give.warnings=FALSE, na.rm=TRUE)$gcv.grid[,1]
# see the 2 column of $gcv.grid to get the effective degress of freedom.
contour( theta,log(lam) , PLL)
## End(Not run)
#plot the candidate matrix by: plot.spatial.design()
plot.spatial.design()
Error in pairs(x$design, ...): error in evaluating the argument 'x' in selecting a method for function 'pairs': Error: argument "x" is missing, with no default
#parallel computing on clusters
library(parallel)
library(Rmpi) #using MPI parallelization scheme
Error in library(Rmpi): there is no package called 'Rmpi'
#setting options
options("scipen" = 10)
parallel::detectCores()
[1] 64
options(mc.cores = parallel::detectCores()) #so that stan and rsampling can use all core
#import data and rasters and plot them
sl_data<- read.csv(file="Sl_AgeCorrectedEx10_300615_2013.csv",header=T,sep=",")
sl_data<-sl_data[,c("ID","Long","Lat","YY","LoAge","UpAge","Urbanizati","tsi","prec","evi","Ex","Pf")]
sl_locs<- read.csv(file="sl_coordinates.csv",header=T,sep=",")
sl_borders<- read.csv(file="boundary_final.csv",header=T,sep=",")
#plot borders and coordinates
coordinates(sl_locs)<-~ Long+Lat
coordinates(sl_borders)<-~ Long+Lat
par(mar=c(0.1,0.1,0.1,0.1))
plot(sl_locs,pch=16,cex=0.7,col="purple")
lines(sl_borders@coords[,"Long"],y= sl_borders@coords[,"Lat"],col= "orange",lwd=2)
#epsg proj4string for sierr-leone
proj_list<-rgdal::make_EPSG()
Error in loadNamespace(name): there is no package called 'rgdal'
proj_list[grep("Sierra",proj_list$note),]
Error in eval(expr, envir, enclos): object 'proj_list' not found
urbanization<-maptools::readAsciiGrid(fname="sl_urb.asc",proj4string=CRS("+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"))
precipitation<-maptools::readAsciiGrid(fname="sl_prec.asc",proj4string=CRS("+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"))
tsi<-maptools::readAsciiGrid(fname="sl_tsi.asc",proj4string=CRS("+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"))
mean_PfPR<-maptools::readAsciiGrid(fname="sle_Correct.asc",proj4string=CRS("+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"))
lapply(list(urbanization,precipitation,mean_PfPR),class)
[[1]]
[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"
[[2]]
[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"
[[3]]
[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"
#checking if the raster properties(projection, dimension,cellsize,extent,ncell) of the rasters agree
urbanization<- as(urbanization,"RasterLayer"); precipitation<- as(precipitation,"RasterLayer");mean_PfPR<- as(mean_PfPR,"RasterLayer"); #changing from "SpatialGridDataFrame"
lapply(list(urbanization,precipitation,mean_PfPR),FUN= proj4string)
[[1]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[2]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[3]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
lapply(list(urbanization,precipitation,mean_PfPR),FUN= proj4string)
[[1]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[2]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[3]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
lapply(list(urbanization,precipitation,mean_PfPR),FUN= proj4string)
[[1]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[2]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
[[3]]
[1] "+init=epsg:2162 +proj=longlat +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84"
#get statistics values from our raster
lapply(list(urbanization,precipitation,mean_PfPR),function(x){cellStats(x, min)})
[[1]]
[1] 0
[[2]]
[1] 72
[[3]]
[1] 0.1299268
lapply(list(urbanization,precipitation,mean_PfPR),function(x){cellStats(x, max)})
[[1]]
[1] 1
[[2]]
[1] 318
[[3]]
[1] 0.7863413
lapply(list(urbanization,precipitation,mean_PfPR),function(x)cellStats(x, range))
[[1]]
[1] 0 1
[[2]]
[1] 72 318
[[3]]
[1] 0.1299268 0.7863413
lapply(list(urbanization,precipitation,mean_PfPR),function(x)quantile(x, probs = c(0.25, 0.75), type=7,names = FALSE))
[[1]]
[1] 0 0
[[2]]
[1] 200 234
[[3]]
[1] 0.4314154 0.5837793
lapply(list(urbanization,precipitation,mean_PfPR),function(x)density(x))
[[1]]
Call:
density.default(x = stats::na.omit(d))
Data: stats::na.omit(d) (63400 obs.); Bandwidth 'bw' = 0.01329
x y
Min. :-0.03987 Min. : 0.0000
1st Qu.: 0.23006 1st Qu.: 0.0000
Median : 0.50000 Median : 0.0000
Mean : 0.50000 Mean : 0.9243
3rd Qu.: 0.76994 3rd Qu.: 0.0000
Max. : 1.03987 Max. :29.3482
[[2]]
Call:
density.default(x = stats::na.omit(d))
Data: stats::na.omit(d) (63400 obs.); Bandwidth 'bw' = 2.501
x y
Min. : 64.5 Min. :0.0000000
1st Qu.:129.7 1st Qu.:0.0000161
Median :195.0 Median :0.0015272
Mean :195.0 Mean :0.0038275
3rd Qu.:260.3 3rd Qu.:0.0049549
Max. :325.5 Max. :0.0214527
[[3]]
Call:
density.default(x = stats::na.omit(d))
Data: stats::na.omit(d) (63696 obs.); Bandwidth 'bw' = 0.01032
x y
Min. :0.09897 Min. :0.000008
1st Qu.:0.27855 1st Qu.:0.069511
Median :0.45813 Median :1.087033
Mean :0.45813 Mean :1.390775
3rd Qu.:0.63771 3rd Qu.:2.580662
Max. :0.81729 Max. :3.586867
freq(urbanization, value=0.09897)
[1] 0
hist(urbanization, main="Distribution of urbanization values", col= "purple", maxpixels=21752940)
hist(precipitation, main="Distribution of precipitation values", col= "purple", maxpixels=21752940)
#A RasterStack is a collection of RasterLayer objects with the same spatial extent and resolution
#raster::extract --ectract values from raster cells
#visualizing our rasters: plot and levelplot is problematic
col <- terrain.colors(5)
image(urbanization, main="Urbanization in SL", col=col)
image(precipitation, main="Precipitation in SL", col=col)
image(mean_PfPR, main="PfPR in SL", col=col)
#Randomly generate points on the sl map
prec<- as(precipitation,"SpatialPolygons")
prec1<-sp::spTransform(x=prec,CRS("+init=epsg:2162 +proj=utm +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84 +units=m +no_defs"))
Error in sp::spTransform(x = prec, CRS("+init=epsg:2162 +proj=utm +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +datum=WGS84 +units=m +no_defs")): package rgdal is required for spTransform methods
(x<-spsample(prec1,n=1000,type="random"))
Error in spsample(prec1, n = 1000, type = "random"): error in evaluating the argument 'x' in selecting a method for function 'spsample': Error: object 'prec1' not found
pts<-sp::spTransform(x=x,CRS(proj4string(prec)))
Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'spTransform' for signature '"matrix", "CRS"'
points(pts,col="red",pch=19,cex=0.4)
Error in points(pts, col = "red", pch = 19, cex = 0.4): object 'pts' not found
#Get Long Lat ranges explicitly from the extent of one of the rasters
extent(urbanization)
class : Extent
xmin : -13.30763
xmax : -10.28263
ymin : 6.928689
ymax : 10.00369
ncell(urbanization)
[1] 133947
xy_new<- xyFromCell(precipitation,1:ncell(precipitation))
colnames(xy_new)<-c("Long","Lat");
#Extract values from raster.
tsi_test<- raster("sl_tsi.asc")
cellFromXY(object=tsi_test,xy=matrix(c(-12.223634,8.899611),ncol=2))
[1] 48047
v<-extract(tsi_test,c(48047),method="simple")
Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'extract' for signature '"RasterLayer"'
xyFromCell(object=tsi_test,cell=48047)
x y
[1,] -12.22013 8.899522
#nd * cd
std_units_design<-scale(x= sl_locs@coords)
sl_locs_center<-attr(std_units_design,"scaled:center")
sl_locs_scale<-attr(std_units_design,"scaled:scale")
knots<- cover.design(R= sl_locs@coords,
nd= 200,#fixed= sl_locs@coords,
nruns= 10,nn=T,num.nn=100,
scale.type="user",R.center=sl_locs_center,R.scale=sl_locs_scale,
P= -Inf,#P= -Inf, Q= Inf gives Minimax design, when P=1,Q=1 simple averaging of distance is used
Q=Inf
#start=sl_locs@coords #we can use "random starting design or sl_locs[c(1,2,3,4,5,65,7,84,13,5,623),]
#DIST= rdist # great circle distance or euclidean
)
Warning in cover.design(R = sl_locs@coords, nd = 200, nruns = 10, nn = T, :
Number of nearst neighbors (nn) reduced to the actual number of candidates
[1] "Greater than 1 optimal design; keeping first one......"
knots_unscaled<- cover.design(R= sl_locs@coords,
nd= 200,#fixed= sl_locs@coords,
nruns= 10,nn=T,num.nn=100,
scale.type="unscaled",R.center=sl_locs_center,R.scale=sl_locs_scale,
P= -Inf,#P= -Inf, Q= Inf gives Minimax design, when P=1,Q=1 simple averaging of distance is used
Q=Inf
#start=sl_locs@coords #we can use "random starting design or sl_locs[c(1,2,3,4,5,65,7,84,13,5,623),]
#DIST= rdist # great circle distance or euclidean
)
Warning in cover.design(R = sl_locs@coords, nd = 200, nruns = 10, nn = T, :
Number of nearst neighbors (nn) reduced to the actual number of candidates
[1] "Greater than 1 optimal design; keeping first one......"
par(mar=c(0.1,0.1,0.1,0.1))
plot(sl_locs,pch=16,cex=0.7,col="purple") ;par(new=T)
lines(sl_borders@coords[,"Long"],y= sl_borders@coords[,"Lat"],col= "orange",lwd=2); par(new=T)
par(new=T)
plot(knots$design,pch="+",col="blue");par(new=T)
plot(knots_unscaled$design,pch="o",col="orange");par(new=T)
