平滑ggplot2中的geom_tile地图 - 插值数据 [英] Smoothing geom_tile map in ggplot2 - Interpolating data

问题描述

我试图从lon-lat数据集绘制空间地图。遵循

这是使用的代码（大部分来自前面引用的问题）

  library（rgdal）
 library（ggplot）
＃读取省份限制
 CV = readOGR（dsn = 。，layer =poligonos_provincia_etrs89）
 CV.pr =子集（CV，CV $ CODINE ==46| CV $ CODINE ==12| CV $ CODINE ==03）
＃将对象转换为data.frame 
 CV2 < -  fortify（CV.pr）
 
＃读取空间数据
 datos.uvi = read.csv（salida -mapa-036.dat，sep =，，header = T，na.strings = -  99.9）
 uvi.temp <-datos.uvi [，c（longitudinal，latitud ，RGlobal）] 
 colnames（uvi.temp）< -c（long，lat，RGlobal）
 
＃plot map 
 ggplot（ ）+ 
 geom_tile（data = uvi.temp，aes（x = long，y = lat，z = RGlobal，fill = RGlobal），alpha = 0.4）+ 
 stat_contour（data = uvi.temp，aes（x = long， y = lat，z = RGlobal））+ 
 geom_path（data = CV2，aes（x = long，y = lat，group = group），color =black，size = 0.6）+ 
 ggtitle（PrevisiónUVI-CV）+ xlab（Longitud）+ ylab（Latitud）+ 
 scale_fill_continuous（name =UVI，low =white，high =red）+ 
 theme_bw（）+ 
 coord_map（）
  

关键是我不' t想要一个瓷砖地图，但一个平滑的梯度地图。运行代码时出现此警告：

警告消息：计算在 stat_contour（）：（list）object
不能被强制输入'double'


我假设这是因为数据请输入下面显示的数据框（带有lon，lat和data的三列）



  dput（uvi.temp）
结构（列表（长= c（-1.25，-1，-1.25，-1，-0.75，-0.5，-0.25，
0,0.25,0.5,0.75，-1.25，-1.25，-1 ，-0.75，-0.5，-0.25,0，
 0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5，
 0.75， - 1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5,0.75，-1.25，
 -1，-0.75，-0.5，-0.25,0,0.25,0.5,0.75， - 1.25，-1，-0.75，
 -0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25，
 0,0.25,0.5 ，0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5，
 0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5 ，0.75，-1.25，
 -1，-0.75，-0.5，-0.25,0,0.25,0.5,0。 75，-1.25，-1，-0.75，
 -0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25，
 0， 0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5，
 0.75，-1.25，-1，-0.75，-0.5，-0.25,0， 0.25,0.5,0.75，-1.25，
 -1，-0.75，-0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75，
 -0.5， - 0.25,0.25,0.5,0.75,0.25,0.1,0.75,0.5,0.25,0 
 0,0.25,0.5,0.75，-1.25，-1，-0.75，-0.5， - 0.25,0,0.25,0.5，
 0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5,0.75，-1.25，
 -1，-0.75， -0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75，
 -0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75， -0.5，-0.25，
 0,0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25,0,0.25,0.5，
 0.75，-1.25，-1 ，-0.75，-0.5，-0.25,0,0.25,0.5,0.75，-1.25，
 -1，-0.75，-0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1 ，-0.75，
 -0.5，-0.25,0,0.25,0.5,0.75，-1.25，-1，-0.75，-0.5，-0.25，
 0,0.25,0.5,0.75）， lat = c（37.75,37.75,37.75,37.75,37.75， 
 37.75,37.75,37.75,37.75,37.75,37.75,38,37.75,37.75,37.75，
 37.75,37.75,37.75,37.75,37.75,37.75,38,38,38,38,38， 
 38,38,38,38,38.25,38.25,38.25,38.25,38.25,38.25,38.25，$ b $ 38.25，38.25，38.5，38.5，38.5，38.5，38.5，38.5，38.5，38.5， 
 38.5，38.75，38.75，38.75，38.75，38.75，38.75，38.75，38.75，
 38.75，39，39，39，39，39，39，39，39，39，39.25，39.25， 39.25，39.25，39.25，39.25，39.25，39.25，39.5，39.5，39.5，39.5，
 39.5，39.5，39.5，39.5，39.5，39.75，39.75，39.75，39.75，39.75， 
 39.75，39.75，39.75，39.75，40，40，40，40，40，40，40，40，40，
 40.25,40.25,40.25,40.25,40.25,40.25,40.25,40.25， 40.25，
 40.5，40.5，40.5，40.5，40.5，40.5，40.5，40.5，40.5，40.75，
 40.75，40.75，40.75，40.75，40.75，40.75，40.75，40.75，37.75，
 37.75,37.75,37.75,37.75,37.75,37.75,37.75,37.75,38,38，$ b $ 38,38,38,38,38,38,38,38.25,38.25,38.25,38.25,38.25， 
 38.25，38.25，38.25， 38.25,38.5,38.5,38.5,38.5,38.5,38.5，
 38.5,38.5,38.5,38.75,38.75,38.75,38.75,38.75,38.75,38.75，$ b $ 38.75,38.75,39,39， 39,39,39,39,39,39,39,39.25,39.25，
 39.25，39.25，39.25，39.25，39.25，39.25，39.25，39.5，39.5，
 39.5，39.5，39.5， 39.5，39.5，39.5，39.5，39.75，39.75，39.75，
 39.75，39.75，39.75，39.75，39.75，39.75，40，40，40，40，40， 40.25,40.25,40.25,40.25,40.25,40.25,40.25，
 40.25,40.25,40.5,40.5,40.5,40.5,40.5,40.5,40.5,40.5，
 40.5,40.75,40.75， 40.75,40.75,40.75,40.75,40.75,40.75，
 40.75），RGlobal = c（469.4264,467.3801,469.4264,467.3801,522.7478，
 443.958,441.283,509.0632,435.3009,432.4717,566.5648,527.1168 ，
 469.4264,467.3801,522.7478,443.958,441.283,509.0632,435.3009，
 432.4717,566.5648,527.1168,459.3489,457.781,439.8998,437.827，
 435.4037,433.0625,430.2155,564.3985,592.1616 ，455.9498,455.622，
 443.1096,43.0866,500美元，欧元，美元，欧元，美元，欧元，美元，欧元，美元，欧元，美元，美元，欧元，美元，欧元，美元，欧元，美元， b 425.5438,419.4045,484.3133,441.4613,445.6218,575.9768,509.4317，
 437.0949,411.3335,412.9341,412.7781,411.7847,432.1784,435.553，
 566.9114,563.8876,414.4346,404.6672,405.9609,405.5721,470.3917 ，
 426.5352，427.227，552.512，584.1269，399.0303，398.1212，398.3982，
 398.6394，462.0748，413.0668，417.8466，481.0621，420.5388，407.8954，
 518.209，390.5623，390.4963，453.643，405.8642 ，405.2392,464.2149，
 413.4183,409.0698,525.5512,447.1689,383.9111,384.0972,400.2987，
 403.8821,405.2475,405.1649,402.3291,519.2192,517.2463,379.0551，
 379.1682，393.9954,396.8607 ，454.2946,394.8534,386.8805,498.197，
 538.8216,387.8813,375.4716,374.1738,377.1193,437.4277,385.8378，
 381.4342,422.0693,365.6135,373.4164,497 .8519,469.4264,467.3801，
 522.7478,443.958,441.283,509.0632,435.3009,432.4717,566.5648，
 527.1168,459.3489,457.781,439.8998,437.827,435.4037,433.0625，
 430.2155,564.3985 ，592.4016,455.9498,455.622,443.1096,43.0866，
 500.2498,428.5327,427.8865,561.918,626.1251,456.2038,458.2242，
 459.7855,431.4439,492.9174,425.223,424.3311,490.1372,450.9101，
 449.7869,586.2458,455.2244,452.9787,514.1946,425.5438,419.4045 
 484.3133,441.4613,445.6218,575.9768,509.4317,437.0949,411.3335 $ b 412.9341,412.7781,411.7847,432.1784,435.553,566.9114,563.8876 ，
 414.4346,404.6672,405.9609,405.5721,470.3917,426.5352,427.227，
 552.512,584.1269,399.0303,398.1212,398.3982,398.6394,462.0748，$ b $ 413.0668,417.8466,481.0621,420.5388,407.8954 ，518.209,390.5623，
 390.4963,453.643,405.8642,405.2392,464.2149,413.4183,409.0698，
 525.5512,447.1689,383.9111,384.0972,400.2987,403.8821,405.247 5，
 405.1649,402.3291,519.2192,517.2463,379.0551,379.1682,393.9954，
 396.8607,454.2946,394.8534,386.8805,498.197,538.8216,387.8813，
 375.4716,374.1738,377.1193,437.4277， 385.8378,381.4342,422.0693，
 365.6135,373.4164,497.8519）），.Names = c（long，lat，RGlobal
），class =data.frame，行。姓名= c（NA，-246L））
  

任何想法？



在此先感谢

解决方案

最后，我可以找到一个解决方案，而不是一个完美的解决方案，一。我会继续寻找更好的插值方法。从

I'm trying to plot an spatial map from a lon-lat data set. Following the example in https://stackoverflow.com/a/19339663/709777 I am able to almost get the desired result

This is the code used (most of it from the previously cited question)

library(rgdal)
library(ggplot)
# read province limits
CV = readOGR(dsn=".", layer="poligonos_provincia_etrs89")
CV.pr=subset(CV, CV$CODINE == "46" | CV$CODINE == "12" | CV$CODINE == "03" )
# convert object to data.frame
CV2 <- fortify(CV.pr)

# read spatial data
datos.uvi=read.csv("salida-mapa-036.dat",sep=",",header=T,na.strings="-99.9")
uvi.temp<-datos.uvi[,c("longitud","latitud","RGlobal")]
colnames(uvi.temp)<-c("long","lat","RGlobal")

# plot map
ggplot() +
  geom_tile(data = uvi.temp, aes(x = long, y = lat, z = RGlobal, fill = RGlobal), alpha = 0.4) +
  stat_contour(data = uvi.temp, aes(x = long, y = lat, z = RGlobal)) +
  geom_path(data = CV2, aes(x = long, y = lat, group = group),color="black", size=0.6) +
  ggtitle("Previsión UVI - CV") +  xlab("Longitud") +  ylab("Latitud") +
  scale_fill_continuous(name = "UVI", low = "white", high = "red") +
  theme_bw() +
  coord_map()

The point is that I don't want a tile map but a smooth gradient map. This warning appears when running the code:

Warning message: Computation failed in stat_contour(): (list) object cannot be coerced to type 'double'

I assume this is because of the data type in the data frame (three columns with lon,lat and data) shown below

dput(uvi.temp)
structure(list(long = c(-1.25, -1, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1.25, -1, -0.75, -0.5, -0.25, 0, 
0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75), lat = c(37.75, 37.75, 37.75, 37.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 37.75, 37.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 38, 38, 38, 38, 
38, 38, 38, 38, 38.25, 38.25, 38.25, 38.25, 38.25, 38.25, 38.25, 
38.25, 38.25, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 
38.5, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 
38.75, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39.25, 39.25, 39.25, 
39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.5, 39.5, 39.5, 39.5, 
39.5, 39.5, 39.5, 39.5, 39.5, 39.75, 39.75, 39.75, 39.75, 39.75, 
39.75, 39.75, 39.75, 39.75, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 
40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.75, 
40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 38, 
38, 38, 38, 38, 38, 38, 38, 38.25, 38.25, 38.25, 38.25, 38.25, 
38.25, 38.25, 38.25, 38.25, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 
38.5, 38.5, 38.5, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 
38.75, 38.75, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39.25, 39.25, 
39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.5, 39.5, 
39.5, 39.5, 39.5, 39.5, 39.5, 39.5, 39.5, 39.75, 39.75, 39.75, 
39.75, 39.75, 39.75, 39.75, 39.75, 39.75, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 
40.25, 40.25, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 
40.5, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 
40.75), RGlobal = c(469.4264, 467.3801, 469.4264, 467.3801, 522.7478, 
443.958, 441.283, 509.0632, 435.3009, 432.4717, 566.5648, 527.1168, 
469.4264, 467.3801, 522.7478, 443.958, 441.283, 509.0632, 435.3009, 
432.4717, 566.5648, 527.1168, 459.3489, 457.781, 439.8998, 437.827, 
435.4037, 433.0625, 430.2155, 564.3985, 592.4016, 455.9498, 455.622, 
443.1096, 431.0866, 500.2498, 428.5327, 427.8865, 561.918, 626.1251, 
456.2038, 458.2242, 459.7855, 431.4439, 492.9174, 425.223, 424.3311, 
490.1372, 450.9101, 449.7869, 586.2458, 455.2244, 452.9787, 514.1946, 
425.5438, 419.4045, 484.3133, 441.4613, 445.6218, 575.9768, 509.4317, 
437.0949, 411.3335, 412.9341, 412.7781, 411.7847, 432.1784, 435.553, 
566.9114, 563.8876, 414.4346, 404.6672, 405.9609, 405.5721, 470.3917, 
426.5352, 427.227, 552.512, 584.1269, 399.0303, 398.1212, 398.3982, 
398.6394, 462.0748, 413.0668, 417.8466, 481.0621, 420.5388, 407.8954, 
518.209, 390.5623, 390.4963, 453.643, 405.8642, 405.2392, 464.2149, 
413.4183, 409.0698, 525.5512, 447.1689, 383.9111, 384.0972, 400.2987, 
403.8821, 405.2475, 405.1649, 402.3291, 519.2192, 517.2463, 379.0551, 
379.1682, 393.9954, 396.8607, 454.2946, 394.8534, 386.8805, 498.197, 
538.8216, 387.8813, 375.4716, 374.1738, 377.1193, 437.4277, 385.8378, 
381.4342, 422.0693, 365.6135, 373.4164, 497.8519, 469.4264, 467.3801, 
522.7478, 443.958, 441.283, 509.0632, 435.3009, 432.4717, 566.5648, 
527.1168, 459.3489, 457.781, 439.8998, 437.827, 435.4037, 433.0625, 
430.2155, 564.3985, 592.4016, 455.9498, 455.622, 443.1096, 431.0866, 
500.2498, 428.5327, 427.8865, 561.918, 626.1251, 456.2038, 458.2242, 
459.7855, 431.4439, 492.9174, 425.223, 424.3311, 490.1372, 450.9101, 
449.7869, 586.2458, 455.2244, 452.9787, 514.1946, 425.5438, 419.4045, 
484.3133, 441.4613, 445.6218, 575.9768, 509.4317, 437.0949, 411.3335, 
412.9341, 412.7781, 411.7847, 432.1784, 435.553, 566.9114, 563.8876, 
414.4346, 404.6672, 405.9609, 405.5721, 470.3917, 426.5352, 427.227, 
552.512, 584.1269, 399.0303, 398.1212, 398.3982, 398.6394, 462.0748, 
413.0668, 417.8466, 481.0621, 420.5388, 407.8954, 518.209, 390.5623, 
390.4963, 453.643, 405.8642, 405.2392, 464.2149, 413.4183, 409.0698, 
525.5512, 447.1689, 383.9111, 384.0972, 400.2987, 403.8821, 405.2475, 
405.1649, 402.3291, 519.2192, 517.2463, 379.0551, 379.1682, 393.9954, 
396.8607, 454.2946, 394.8534, 386.8805, 498.197, 538.8216, 387.8813, 
375.4716, 374.1738, 377.1193, 437.4277, 385.8378, 381.4342, 422.0693, 
365.6135, 373.4164, 497.8519)), .Names = c("long", "lat", "RGlobal"
), class = "data.frame", row.names = c(NA, -246L))

Any idea?

Thanks in advance

解决方案

Finally I could find a solution, not a perfect one but a good enough one. I'll keep searching for better interpolation methods. From the question in https://gis.stackexchange.com/q/169184/9227 I made this code

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
library(rgdal)

# Reading three data frames
datos.uvi.1 <- read.csv(file = "./salida-mapa-012.dat",header = TRUE)
datos.uvi.2 <- read.csv(file = "./salida-mapa-036.dat",header = TRUE)
datos.uvi.3 <- read.csv(file = "./salida-mapa-060.dat",header = TRUE)

# Looking for RGlobal max
new.uvi <- data.frame(datos.uvi.1$RGlobal,datos.uvi.2$RGlobal,datos.uvi.3$RGlobal)
max=apply(new.uvi, 1, max, na.rm=FALSE)
datos.uvi=data.frame(datos.uvi.1$longitud,datos.uvi.1$latitud,max)
colnames(datos.uvi)<-c("longitud","latitud","RGlobal")

# Define x & y as longitude and latitude
datos.uvi$x <- datos.uvi$longitud 
datos.uvi$y <- datos.uvi$latitud

coordinates(datos.uvi) = ~x + y

# Reading shapefiles for province limits (black and blue limits in the map)
provincias = readOGR(dsn=".", layer="poligonos_provincia_etrs89")
pr1=subset(provincias, provincias$CODINE == "46" | provincias$CODINE == "12" | provincias$CODINE == "03" )
pr2=subset(provincias, provincias$CODINE == "43" | provincias$CODINE == "16" | provincias$CODINE == "30"  | provincias$CODINE == "02" | provincias$CODINE == "50" )
pr1 <- fortify(pr1)
pr2 <- fortify(pr2)

# Interpolation area
x.range <- as.numeric(c(-1.75, 1))  # min/max longitude
y.range <- as.numeric(c(37.5, 41))  # min/max latitude

# Create a gridded structure
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], to = y.range[2], by = 0.1))
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

#Interpolate surface and fix the output. Apply idw model for the data
idw <- idw(formula = RGlobal ~ 1, locations = datos.uvi, newdata = grd)  
idw.output = as.data.frame(idw) 
names(idw.output)[1:3] <- c("long", "lat", "RGlobalmax")

# Plot
  ggplot() + geom_tile(data = idw.output, alpha = 0.8, aes(x = long, y = lat, fill = RGlobalmax)) +
    scale_fill_gradient(low = "cyan", high = "orange",name = "UVI") +
    geom_path(data = pr2, aes(long, lat, group = group), colour = "blue") +
    geom_path(data = pr1, aes(long, lat, group = group), colour = "black") +
    coord_map(xlim = c(-1.7, 1),ylim = c(37.6,40.9)) +
    ggtitle("Previsión UVI - DD/MM/YYYY") + xlab(" ") + ylab(" ")

And the final output

这篇关于平滑ggplot2中的geom_tile地图 - 插值数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！