r - 将火灾数据作为具有纬度和经度值的点添加到使用 R 中的 OpenAir 包生成的气团回弹轨迹图中
问题描述
我目前有一系列气团后向轨迹图,其中按后向轨迹交叉的频率网格,使用 R 的 OpenAir 包生成。我想在这些图中添加指示火灾位置的点。从 .csv 文件中读取的点,其中每个独特的火灾都有一个纬度和经度值。
我已阅读 OpenAir 手册,但没有看到任何解决方案。我尝试使用以下代码添加要点:
Week_One_Freq + points (FireData$longitude, FireData$latitude, col = "red", cex = .6)
但是遇到错误:
Error in Week_One_Freq + points(FireData$longitude, FireData$latitude,: non-numeric argument to binary operator
我已经测试了问题是否与火灾数据有关。我绘制的点如下:
library(rworldmap)
newmap <- getMap(resolution = "low")
plot(newmap, xlim = c(-80, 80), ylim = c(-40, 40), asp = 1)
points (FireData$longitude, FireData$latitude, col = "red", cex = .6)
它工作正常,正确绘制火灾位置。
这是我的频率轨迹图代码:
Week_One_Freq <- trajLevel(Week_1_Traj, statistic = "frequency",
smooth = TRUE,
map.cols = "grey40", map.alpha = 0.2,
projection = "mercator", parameters = NULL,
grid.col = "transparent",
key.pos = "right",
orientation = c(90,0,0),
npoints = 20, origin = TRUE,
xlim= c(-80, 80), ylim = c(-40, 40))
到目前为止,我可以在地图上绘制火灾数据点,并且可以为气团返回轨迹生成一系列频率图。我想将两者结合起来,以便在一张地图上同时显示返回轨迹和指示火灾位置的点。
编辑:用于创建频率图的轨迹数据示例:
structure(list(X = 1:5, Receptor = c(1L, 1L, 1L, 1L,
1L, 1L), Year = c(2017L, 2017L, 2017L, 2017L, 2017L, 2017L),
Month = c(10L, 10L, 10L, 10L, 10L, 10L), Day = c(6L, 5L,
5L, 5L, 5L, 5L), Hour = c(0L, 23L, 22L, 21L, 20L, 19L), hour.inc = 0:-5,
lat = c(-3.162, -3.189, -2.937, -2.87, -2.809, -2.753), lon = c(-60.6,
-60.439, -60.282, -60.13, -60.986, -60.851), height = c(500,
502.2, 500.6, 495.2, 486.4, 474.3), pressure = c(944.2, 942.9,
942.4, 942.3, 942.7, 943.5), date = structure(c(1L, 1L, 1L,
1L, 1L, 1L), .Label = c("06/10/2017 00:00", "06/10/2017 03:00",
"06/10/2017 06:00", "06/10/2017 09:00", "06/10/2017 12:00",
"06/10/2017 15:00", "06/10/2017 18:00", "06/10/2017 21:00",
"07/10/2017 00:00", "07/10/2017 03:00", "07/10/2017 06:00",
"07/10/2017 09:00", "07/10/2017 12:00", "07/10/2017 15:00",
"07/10/2017 18:00", "07/10/2017 21:00", "08/10/2017 00:00",
"08/10/2017 03:00", "08/10/2017 06:00", "08/10/2017 09:00",
"08/10/2017 12:00", "08/10/2017 15:00", "08/10/2017 18:00",
"08/10/2017 21:00", "09/10/2017 00:00", "09/10/2017 03:00",
"09/10/2017 06:00", "09/10/2017 09:00", "09/10/2017 12:00",
"09/10/2017 15:00", "09/10/2017 18:00", "09/10/2017 21:00",
"10/10/2017 00:00", "10/10/2017 03:00", "10/10/2017 06:00",
"10/10/2017 09:00", "10/10/2017 12:00", "10/10/2017 15:00",
"10/10/2017 18:00", "10/10/2017 21:00", "11/10/2017 00:00",
"11/10/2017 03:00", "11/10/2017 06:00", "11/10/2017 09:00",
"11/10/2017 12:00", "11/10/2017 15:00", "11/10/2017 18:00",
"11/10/2017 21:00", "12/10/2017 00:00", "12/10/2017 03:00",
"12/10/2017 06:00", "12/10/2017 09:00", "12/10/2017 12:00",
"12/10/2017 15:00", "12/10/2017 18:00", "12/10/2017 21:00"
), class = "factor"), Date = structure(c(1L, 1L, 1L, 1L,
1L, 1L), .Label = c("06/10/2017", "07/10/2017", "08/10/2017",
"09/10/2017", "10/10/2017", "11/10/2017", "12/10/2017"), class =
"factor")), row.names = c(NA,
6L), class = "data.frame")
火灾数据示例:
structure(list(acq_date = structure(c(6L, 6L, 6L, 6L, 6L, 6L), .Label =
c("01/11/2017", "02/11/2017", "03/11/2017", "04/11/2017", "05/11/2017",
"06/10/2017", "06/11/2017", "07/10/2017", "08/10/2017", "09/10/2017",
"10/10/2017", "11/10/2017", "12/10/2017", "13/10/2017", "14/10/2017",
"15/10/2017", "16/10/2017", "17/10/2017", "18/10/2017", "19/10/2017",
"20/10/2017", "21/10/2017", "22/10/2017", "23/10/2017", "24/10/2017",
"25/10/2017", "26/10/2017", "27/10/2017", "28/10/2017", "29/10/2017",
"30/10/2017", "31/10/2017"), class = "factor"), latitude = c(-0.18499,
-0.18777, -0.18832, -0.1932, -0.19654, -0.21077), longitude = c(-48.89484,
-48.89885, -48.89536, -49.06121, -49.06173, -48.97051)), row.names = c(NA,
6L), class = "data.frame")
解决方案
有人在https://github.com/davidcarslaw/openair/issues/68解决了同样的问题。您可以使用 latticeExtra 中的图层函数将点添加到最后一个格子对象,但首先您需要使用 mapproject 将您的点放到正确的投影上。我已经扩展了他们的示例以包含一些附加信息,例如所需的库。
library(openair)
library(maps)
library(mapproj)
library(latticeExtra)
# import your data
Trajectory_Data <- read.table("Analysis.txt")
Fire_Data <- read.csv("Fire_Data.csv")
lon <- Fire_Data$longitude
lat <- Fire_Data$latitude
#convert points to projected lat/long
tmp <- mapproject(x=lon, y=lat, projection='lambert', parameters = c(51, 51), orientation = c(90, 0, 0))
# I found the projection, parameters, and orientation by looking at the default values for arguments in the trajLevel function
# Plot trajectories
trajLevel(subset(Trajectory_Data, lon > -70 & lon < -20 & lat > -15 & lat < 15),
pollutant = "BC", statistic = "cwt", smooth = TRUE, col = "increment", border = NA)
# add the points. you can adjust the size of the points with cex
c <- trellis.last.object() + layer(lpoints(tmp$x, tmp$y, pch=16, cex=0.05, col="black"))
print(c)
推荐阅读
- amazon-cloudformation - Cloudformation aws-glue 内联命令
- python - Instagram API:查找特定半径内所有图像的纬度/经度坐标
- c - 长字符串上的 Segfault 11,在访问字符串之前,仅当字符串 > 14 时
- php - 从 mysql 回显特殊字符
- batch-file - Cleanmgr 使用命令或批处理
- c# - 从 HttpPostedFileBase 转换为图像产生错误由于位置
- php - 如何获取 jQuery 对象值?
- pandas - 熊猫数据框和重复名称
- c - 如果有多个 return 语句,如何在 C 中执行 return 语句?
- python - Python递归仅在必要时评估