在 R 中对齐大栅格和空间点
Aligning large raster and SpatialPoints in R
我有两个数据集:
- 一个栅格文件,其中包含澳大利亚的每像素植被分类。
该文件被描述为具有 EPSG:3577 https://epsg.io/3577. All additional meta data, along with data itself can be found at http://www.agriculture.gov.au/abares/forestsaustralia/forest-data-maps-and-tools/spatial-data/forest-cover
的投影
- 一个 csv 文件 包含纬度列、经度列和各种其他变量的数据。
由于这也是一个大数据集,我将只提供它的头部:
firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083,
-23.5908, -23.689), longitude = c(147.6421, 148.8484, 148.1897,
147.2966, 150.1676, 150.0994), brightness = c(320.7, 316.7, 320.5,
320.1, 315.8, 314.3), scan = c(1.6, 2.1, 1.8, 1.7, 2.7, 2.7),
track = c(1.3, 1.4, 1.3, 1.3, 1.6, 1.6), acq_date = c("2017-01-01",
"2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01"
), acq_time = c(47L, 47L, 47L, 47L, 47L, 47L), satellite = c("Terra",
"Terra", "Terra", "Terra", "Terra", "Terra"), instrument = c("MODIS",
"MODIS", "MODIS", "MODIS", "MODIS", "MODIS"), confidence = c(34L,
26L, 36L, 53L, 33L, 22L), version = c(6.2, 6.2, 6.2, 6.2,
6.2, 6.2), bright_t31 = c(299.6, 295.3, 299.7, 296.6, 291.7,
289.3), frp = c(19, 25.8, 22.9, 17.6, 35.4, 30), daynight = c("D",
"D", "D", "D", "D", "D"), type = c(0L, 0L, 0L, 0L, 0L, 0L
)), row.names = c(NA, -6L), class = c("data.table", "data.frame"
), .internal.selfref = <pointer: 0x10200e2e0>)
我的目标很简单:
绘制栅格数据,然后将 csv 文件中的每一行绘制为栅格顶部的一个点。
我对 plot()
、levelplot()
或 ggplot()
很满意,尽管后者看起来很愚蠢,因为将这个巨大的栅格转换为 data.frame 需要太多时间.
到目前为止,我的代码如下所示:
library(rgdal)
library(raster)
library(rasterVis)
library(sp)
- 加载光栅文件
test <- raster('w001000.adf')
- 然后为 data.frame 创建坐标并使其成为 SpatialPoints 格式
coordinates(firms) <- ~longitude + latitude
firms_pts <- SpatialPoints(coords=coordinates(firms), proj4string = CRS("+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
- 将点数据转换为与光栅
相同的crs
firms_lyr <- spTransform(firms_pts, crs(test))
当我尝试绘制此时的数据时,我看不到任何点。
- 所以我在两个数据集中使用相同的范围,然后大致在彼此之上绘制。
extent(test) <- extent(firms)
- 情节
plot(test)
points(firms_lyr, pch=16, cex=0.1)
现在...烦人的是,数据集在空间上没有正确对齐
显然这与投影有关。
我很确定我在步骤 3 中做错了什么 and/or 4
这有效吗?我不确定这些点是不是在远东。
数据
firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083, -23.5908, -23.689),
longitude = c(147.6421, 148.8484, 148.1897, 147.2966, 150.1676, 150.0994)),
.Names = c("longitude", "latitude"),
class = "data.frame", row.names = c(NA, -6L))
xy <- firms[,c(2,1)]
firms_pts <- SpatialPointsDataFrame(coords = xy, data = firms,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
使用 adf 的 Shapefile 重投影
library(raster)
setwd(paste(dirname(rstudioapi::getActiveDocumentContext()$path), "/aus_for18", sep=""))
test <- raster("aus_for18.adf")
firms_lyr <- spTransform(firms_pts, crs(test))
绘图
plot(test)
points(firms_lyr, pch=16, cex=1) #Queensland points
我有两个数据集:
- 一个栅格文件,其中包含澳大利亚的每像素植被分类。
该文件被描述为具有 EPSG:3577 https://epsg.io/3577. All additional meta data, along with data itself can be found at http://www.agriculture.gov.au/abares/forestsaustralia/forest-data-maps-and-tools/spatial-data/forest-cover
的投影- 一个 csv 文件 包含纬度列、经度列和各种其他变量的数据。
由于这也是一个大数据集,我将只提供它的头部:
firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083,
-23.5908, -23.689), longitude = c(147.6421, 148.8484, 148.1897,
147.2966, 150.1676, 150.0994), brightness = c(320.7, 316.7, 320.5,
320.1, 315.8, 314.3), scan = c(1.6, 2.1, 1.8, 1.7, 2.7, 2.7),
track = c(1.3, 1.4, 1.3, 1.3, 1.6, 1.6), acq_date = c("2017-01-01",
"2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01", "2017-01-01"
), acq_time = c(47L, 47L, 47L, 47L, 47L, 47L), satellite = c("Terra",
"Terra", "Terra", "Terra", "Terra", "Terra"), instrument = c("MODIS",
"MODIS", "MODIS", "MODIS", "MODIS", "MODIS"), confidence = c(34L,
26L, 36L, 53L, 33L, 22L), version = c(6.2, 6.2, 6.2, 6.2,
6.2, 6.2), bright_t31 = c(299.6, 295.3, 299.7, 296.6, 291.7,
289.3), frp = c(19, 25.8, 22.9, 17.6, 35.4, 30), daynight = c("D",
"D", "D", "D", "D", "D"), type = c(0L, 0L, 0L, 0L, 0L, 0L
)), row.names = c(NA, -6L), class = c("data.table", "data.frame"
), .internal.selfref = <pointer: 0x10200e2e0>)
我的目标很简单:
绘制栅格数据,然后将 csv 文件中的每一行绘制为栅格顶部的一个点。
我对 plot()
、levelplot()
或 ggplot()
很满意,尽管后者看起来很愚蠢,因为将这个巨大的栅格转换为 data.frame 需要太多时间.
到目前为止,我的代码如下所示:
library(rgdal)
library(raster)
library(rasterVis)
library(sp)
- 加载光栅文件
test <- raster('w001000.adf')
- 然后为 data.frame 创建坐标并使其成为 SpatialPoints 格式
coordinates(firms) <- ~longitude + latitude
firms_pts <- SpatialPoints(coords=coordinates(firms), proj4string = CRS("+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
- 将点数据转换为与光栅 相同的
crs
firms_lyr <- spTransform(firms_pts, crs(test))
当我尝试绘制此时的数据时,我看不到任何点。
- 所以我在两个数据集中使用相同的范围,然后大致在彼此之上绘制。
extent(test) <- extent(firms)
- 情节
plot(test)
points(firms_lyr, pch=16, cex=0.1)
现在...烦人的是,数据集在空间上没有正确对齐
显然这与投影有关。
我很确定我在步骤 3 中做错了什么 and/or 4
这有效吗?我不确定这些点是不是在远东。
数据
firms <- structure(list(latitude = c(-20.5897, -22.4119, -21.132, -23.9083, -23.5908, -23.689),
longitude = c(147.6421, 148.8484, 148.1897, 147.2966, 150.1676, 150.0994)),
.Names = c("longitude", "latitude"),
class = "data.frame", row.names = c(NA, -6L))
xy <- firms[,c(2,1)]
firms_pts <- SpatialPointsDataFrame(coords = xy, data = firms,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
使用 adf 的 Shapefile 重投影
library(raster)
setwd(paste(dirname(rstudioapi::getActiveDocumentContext()$path), "/aus_for18", sep=""))
test <- raster("aus_for18.adf")
firms_lyr <- spTransform(firms_pts, crs(test))
绘图
plot(test)
points(firms_lyr, pch=16, cex=1) #Queensland points