Reputation: 101
I'm trying to do a universal kriging to get a plot with the predicted values inside a map of a data state. However, after fitting the model and calling the predict.gstat
function, the var1.pred
values are all NA
. It is likely that the error is in my grid. However, my theoretical knowledge is minimal to try to understand the reason for these results. I've been studying the subject for two days and looking for examples of universal kriging on the internet. However, there is very little on the subject. The vast majority of literature works with ordinary kriging. Here's an MRC for someone to help me.
library(geobr)
library(sf)
library(automap)
library(gstat)
ms <- read_state(code_state = 'MS')
newdados=structure(list(ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68), lat = c(-19.7584844917999, -19.7554756324444,
-19.8092322281484, -19.7431342667974, -19.7796753287081, -19.897124623193,
-20.2657699257019, -20.2553585407888, -20.3763573198698, -20.4701579483352,
-19.1701011565229, -20.6114105188201, -19.8137313028203, -19.7735844880099,
-19.7260209241669, -19.7395098284587, -19.8578624930131, -19.8127360438838,
-20.4818425644023, -20.4323762329049, -20.5532171100791, -20.6365778239683,
-20.5380236231899, -20.3542398943282, -20.1380771533406, -20.1035964386219,
-20.378798539174, -20.0043199526624, -20.0021734087336, -21.7723132812797,
-21.0631910998205, -20.0381319640795, -20.242115060518, -20.4691114712285,
-20.689567270837, -21.2401609713894, -20.5661396309294, -19.5895515580745,
-19.5324391256981, -19.5204992041817, -19.3777493035806, -19.9814972284517,
-19.8558116452559, -19.9557823498238, -19.955790549824, -19.9903814107959,
-19.8740324100993, -19.9608918567551, -21.0323058296415, -20.8374435298574,
-21.1000506253652, -20.5371714434305, -20.5951455682202, -19.5209019051429,
-20.1936126002268, -20.4491723280869, -20.2987424493016, -20.2992842851162,
-20.3148633451842, -20.4273136654749, -20.242115060518, -20.689567270837,
-18.9634135835693, -18.2301017520185, -17.5313680545638, -20.5951394836648,
-19.7584844002068, -21.5591240076556), long = c(-54.5479795414318,
-54.545879425965, -54.5477176103294, -54.5056103200482, -54.6423494701182,
-54.6609790609671, -54.790257377105, -54.8521503579494, -54.6681786167462,
-54.4731603631335, -54.4328582381536, -54.3241700992432, -54.8460546539279,
-54.8550862580805, -54.8724328753037, -54.8697263419886, -55.1754996795801,
-54.8451614490949, -55.4315749658576, -55.3980560763053, -55.4823038247501,
-55.4098136881897, -55.444735756653, -55.3711721363094, -54.4007870407057,
-54.4280636148514, -54.3712754373255, -54.7765230337685, -54.7781965398355,
-53.750890074785, -54.1242982059117, -53.7984344810548, -53.0707669379434,
-53.3783621739118, -53.3847667416623, -53.2180476520898, -52.8052796761472,
-55.1243857182411, -54.9977852650296, -54.9408303088306, -55.1213025096583,
-54.4561842044892, -54.8528469391583, -54.8660640019876, -54.8637427647631,
-54.7871036213883, -54.8013158405162, -54.8293721843903, -55.2693799558283,
-54.920123671528, -55.0679588584736, -55.0472301343307, -55.3523474341436,
-54.6774751240358, -54.9422256927908, -55.047367149261, -54.9077553865345,
-54.9092606584648, -54.8978976652863, -54.9792504909982, -53.0707669379434,
-53.3847667416623, -54.509658344922, -53.4689512979018, -54.73766298066,
-55.3523440459297, -54.5479794474474, -54.5144340910516), nota = c(58.9026453074775,
43.0734389907771, 56.5105952766095, 36.5680007243752, 43.2, 43.1280350107967,
48.9395514223195, 48.3213464192609, 89.75, 55, 19, 48.0043028169014,
49.622439256789, 45.49609375, 43.2, 43.1710579651398, 37.2179607109448,
50, 48.0003632401017, 45.6276812921986, 42.3627225578784, 48.925,
47.055, 38.6377163904236, 46.5491855979426, 51.8186416184971,
41.9029411764706, 55, 55, 43.2001071319253, 40.4316425648022,
43.1690522819569, 38.88, 44, 39.8602471239881, 41.8588324106694,
41.9094516038661, 49.0645161290323, 58.6666666666667, 43.8720472440945,
40.5, 59.3194635488308, 55, 37.5, 37.5, 55, 48.925, 60, 43.2570422535211,
43.2, 76, 46.35, 46.7787514906126, 53.8175, 49.4354666982136,
48.09375, 45.7303921568627, 43.2, 35.4568965517241, 35.625, 41.76,
44.8427780144866, 15, 28.695652173913, 68.8349514563107, 62.0808951995341,
73.7763345101965, 41.6854924202269), dista = c(30, 31, 23, 32,
37, 32, 33, 36, 8, 8, 40, 55, 3, 9, 16, 15, 15, 2, 31, 40, 31,
22, 28, 51, 1, 1, 51, 14, 13, 85, 105, 45, 35, 40, 65, 52, 20,
28, 11, 11, 32, 20, 4, 2, 3, 12, 6, 6.7, 44, 13, 31, 28, 14,
16, 32, 21, 22, 20, 20, 14, 35, 65, 75, 27, 30, 12.5, 25, 25),
logvalor = c(3.40119738166216, 3.43398720448515, 3.13549421592915,
3.46573590279973, 3.61091791264422, 3.46573590279973, 3.49650756146648,
3.58351893845611, 2.07944154167984, 2.07944154167984, 3.68887945411394,
4.00733318523247, 1.09861228866811, 2.19722457733622, 2.77258872223978,
2.70805020110221, 2.70805020110221, 0.693147180559945, 3.43398720448515,
3.68887945411394, 3.43398720448515, 3.09104245335832, 3.3322045101752,
3.93182563272433, 0, 0, 3.93182563272433, 2.63905732961526,
2.56494935746154, 4.44265125649032, 4.65396035015752, 3.80666248977032,
3.55534806148941, 3.68887945411394, 4.17438726989564, 3.95124371858143,
2.99573227355399, 3.3322045101752, 2.39789527279837, 2.39789527279837,
3.46573590279973, 2.99573227355399, 1.38629436111989, 0.693147180559945,
1.09861228866811, 2.484906649788, 1.79175946922805, 1.90210752639692,
3.78418963391826, 2.56494935746154, 3.43398720448515, 3.3322045101752,
2.63905732961526, 2.77258872223978, 3.46573590279973, 3.04452243772342,
3.09104245335832, 2.99573227355399, 2.99573227355399, 2.63905732961526,
3.55534806148941, 4.17438726989564, 4.31748811353631, 3.29583686600433,
3.40119738166216, 2.52572864430826, 3.2188758248682, 3.2188758248682
)), row.names = c(NA, -68L), class = "data.frame")
dadosop3 <- st_as_sf(newdados,coords = c('long','lat'),crs=4674)
dadspa <- as(dadosop3,'Spatial')
grade.ms <- st_make_grid(ms,cellsize=c(0.2,0.2))%>%
st_as_sf()
grade.ms <- grade.ms[st_contains(ms,grade.ms,sparse=FALSE),]
plot(grade.ms)
grade.ms$nota <- seq(1,100,l=696)
grade.ms$dista <- seq(1,115,l=696)
spgrspa <- as(grade.ms,'Spatial')
vario2 <- autofitVariogram(logvalor ~ nota+dista,dadspa)
modelo <- gstat(formula = logvalor ~nota+dista,
model = vario2$var_model,
data = dadspa)
valor_int <- predict(modelo,spgrspa)%>%
st_as_sf()
Upvotes: 1
Views: 144