将坐标、半径和站点类型数据从 .csv 转换为 R 中的栅格
Converting coordinates, radius and site type data from .csv to raster in R
我在将 .csv 文件转换为 R 中的光栅时遇到了一些问题...我的 .csv 文件包含坐标（经度和纬度）、半径（以度为单位）和站点类型。我能够将坐标转换为光栅并能够使用 st_buffer() 绘制圆，但我面临两个问题：
- 我无法将圆圈转换成光栅...我尝试使用
rasterize() 和 fasterize() 两者都不起作用我得到的只是一个空光栅层
- 我好像不能根据站点类型对坐标和圆圈进行分类
知道我可能做错了什么吗？以及如何对我的圈子进行分类？
提前致谢！
这是我使用的代码：
> head(sp_csv_data)
Longitude Latitude Radius Site_Type
1 -177.87567 -24.715167 10 MIG
2 -83.21360 14.401800 1 OBS
3 -82.59392 9.589192 1 NES
4 -82.41060 9.492750 1 NES;BRE
5 -81.17555 7.196750 5 OBS
6 -80.95770 8.852700 1 NES
##Projection systems used
rob_pacific <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" # Best to define these first so you don't make mistakes below
longlat <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
####Converting to raster####
# Creating a empty raster at 0.5° resolution (you can increase the resolution to get a better border precision)
rs <- raster(ncol = 360*2, nrow = 180*2)
rs[] <- 1:ncell(rs)
crs(rs) <- CRS(longlat)
##Converting to raster
sp_raster <- rasterize(sp_csv_data[,1:2], rs, sp_csv_data[,3])
# Resampling to make sure that it's in the same resolution as sampling area
sp_raster <- resample(sp_raster, rs, resample = "ngb")
#converting into an sf spatial polygon dataframe
sp_raster <- as(sp_raster, "SpatialPolygonsDataFrame")
species_sp <- spTransform(sp_raster, CRS(longlat))
# Define a long & slim polygon that overlaps the meridian line & set its CRS to match that of world
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90),
c(0, 90),
c(0, -90),
c(-0.0001, -90),
c(-0.0001, 90)))) %>%
st_sfc() %>%
st_set_crs(longlat)
# Transform the species distribution polygon object to a Pacific-centred projection polygon object
sp_robinson <- species_sp %>%
st_as_sf() %>%
st_difference(polygon) %>%
st_transform(crs = rob_pacific)
# There is a line in the middle of Antarctica. This is because we have split the map after reprojection. We need to fix this:
bbox1 <- st_bbox(sp_robinson)
bbox1[c(1,3)] <- c(-1e-5,1e-5)
polygon1 <- st_as_sfc(bbox1)
crosses1 <- sp_robinson %>%
st_intersects(polygon1) %>%
sapply(length) %>%
as.logical %>%
which
# Adding buffer 0
sp_robinson[crosses1, ] %<>%
st_buffer(0)
# Adding the circles to the coordinates
sp_robinson2 <- st_buffer(sp_robinson, dist = radius)
> print(sp_robinson2)
Simple feature collection with 143 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -17188220 ymin: -5706207 xmax: 17263210 ymax: 6179000
CRS: +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
First 10 features:
layer geometry
1 5 POLYGON ((3556791 4766657, ...
2 10 POLYGON ((13713529 4995696,...
3 10 POLYGON ((12834403 4946927,...
4 10 POLYGON ((9991443 4801974, ...
5 5 POLYGON ((4254202 4304190, ...
6 5 POLYGON ((11423719 4327354,...
7 10 POLYGON ((9582710 4282247, ...
8 10 POLYGON ((588877.2 4166512,...
9 5 POLYGON ((4522824 3894919, ...
10 10 POLYGON ((3828685 3886205, ...
sp_robinson3 <- fasterize(sp_robinson2, rs)
> print(sp_robinson3)
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
crs : +proj=robin +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : NA, NA (min, max)
我想将 sp_robinson2 转换成一个名为 sp_robinson3 的栅格，但如您所见，fasterize() 和 rasterize() 都给我一个空栅格图层...
最后光栅化不行的原因很明显：矢量和光栅的crs不匹配。但是你能稍微编辑一下你的问题来解释你想要实现的目标吗？创建一个栅格然后创建多边形然后再次栅格化这些是非常奇怪的。我的印象是你让事情变得比需要的复杂得多。你也谈到圈子。哪些圈子？我猜你可能想围绕你的观点画圈，但这不是你在做什么。逐步弄清楚事情可能会有所帮助，首先弄清楚如何获得您想要的一般结果，然后如何使其以太平洋为中心。
以下是代码第一部分的清理版本。它还使它具有可重现性。您需要在代码中创建示例，如下所示：
lon <- c(-177.87567, -83.2136, -82.59392, -82.4106, -81.17555, -80.9577)
lat <- c(-24.715167, 14.4018, 9.589192, 9.49275, 7.19675, 8.8527)
radius <- c(10, 1, 1, 1, 5, 1)
site <- c("MIG", "OBS", "NES", "NES;BRE", "OBS", "NES")
sp_csv_data <- data.frame(lon, lat, radius, site)
## cleaned up projection definitions
rob_pacific <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
longlat <- "+proj=longlat +datum=WGS84"
##Converting to raster
# Creating a empty raster at 0.5° resolution
rs <- raster(res=0.5, crs=lonlat)
values(rs) <- 1:ncell(rs)
sp_raster <- rasterize(sp_csv_data[,1:2], rs, sp_csv_data[,3])
## makes no sense: sp_raster <- resample(sp_raster, rs, resample = "ngb")
#converting into an SpatialPolygonsDataframe
species_sp <- as(sp_raster, "SpatialPolygonsDataFrame")
## makes no sense: species_sp <- spTransform(sp_raster, CRS(longlat))
我在将 .csv 文件转换为 R 中的光栅时遇到了一些问题...我的 .csv 文件包含坐标（经度和纬度）、半径（以度为单位）和站点类型。我能够将坐标转换为光栅并能够使用 st_buffer() 绘制圆，但我面临两个问题：
- 我无法将圆圈转换成光栅...我尝试使用
rasterize()和fasterize()两者都不起作用我得到的只是一个空光栅层 - 我好像不能根据站点类型对坐标和圆圈进行分类
知道我可能做错了什么吗？以及如何对我的圈子进行分类？ 提前致谢！
这是我使用的代码：
> head(sp_csv_data)
Longitude Latitude Radius Site_Type
1 -177.87567 -24.715167 10 MIG
2 -83.21360 14.401800 1 OBS
3 -82.59392 9.589192 1 NES
4 -82.41060 9.492750 1 NES;BRE
5 -81.17555 7.196750 5 OBS
6 -80.95770 8.852700 1 NES
##Projection systems used
rob_pacific <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" # Best to define these first so you don't make mistakes below
longlat <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
####Converting to raster####
# Creating a empty raster at 0.5° resolution (you can increase the resolution to get a better border precision)
rs <- raster(ncol = 360*2, nrow = 180*2)
rs[] <- 1:ncell(rs)
crs(rs) <- CRS(longlat)
##Converting to raster
sp_raster <- rasterize(sp_csv_data[,1:2], rs, sp_csv_data[,3])
# Resampling to make sure that it's in the same resolution as sampling area
sp_raster <- resample(sp_raster, rs, resample = "ngb")
#converting into an sf spatial polygon dataframe
sp_raster <- as(sp_raster, "SpatialPolygonsDataFrame")
species_sp <- spTransform(sp_raster, CRS(longlat))
# Define a long & slim polygon that overlaps the meridian line & set its CRS to match that of world
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90),
c(0, 90),
c(0, -90),
c(-0.0001, -90),
c(-0.0001, 90)))) %>%
st_sfc() %>%
st_set_crs(longlat)
# Transform the species distribution polygon object to a Pacific-centred projection polygon object
sp_robinson <- species_sp %>%
st_as_sf() %>%
st_difference(polygon) %>%
st_transform(crs = rob_pacific)
# There is a line in the middle of Antarctica. This is because we have split the map after reprojection. We need to fix this:
bbox1 <- st_bbox(sp_robinson)
bbox1[c(1,3)] <- c(-1e-5,1e-5)
polygon1 <- st_as_sfc(bbox1)
crosses1 <- sp_robinson %>%
st_intersects(polygon1) %>%
sapply(length) %>%
as.logical %>%
which
# Adding buffer 0
sp_robinson[crosses1, ] %<>%
st_buffer(0)
# Adding the circles to the coordinates
sp_robinson2 <- st_buffer(sp_robinson, dist = radius)
> print(sp_robinson2)
Simple feature collection with 143 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -17188220 ymin: -5706207 xmax: 17263210 ymax: 6179000
CRS: +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
First 10 features:
layer geometry
1 5 POLYGON ((3556791 4766657, ...
2 10 POLYGON ((13713529 4995696,...
3 10 POLYGON ((12834403 4946927,...
4 10 POLYGON ((9991443 4801974, ...
5 5 POLYGON ((4254202 4304190, ...
6 5 POLYGON ((11423719 4327354,...
7 10 POLYGON ((9582710 4282247, ...
8 10 POLYGON ((588877.2 4166512,...
9 5 POLYGON ((4522824 3894919, ...
10 10 POLYGON ((3828685 3886205, ...
sp_robinson3 <- fasterize(sp_robinson2, rs)
> print(sp_robinson3)
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
crs : +proj=robin +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : NA, NA (min, max)
我想将 sp_robinson2 转换成一个名为 sp_robinson3 的栅格，但如您所见，fasterize() 和 rasterize() 都给我一个空栅格图层...
最后光栅化不行的原因很明显：矢量和光栅的crs不匹配。但是你能稍微编辑一下你的问题来解释你想要实现的目标吗？创建一个栅格然后创建多边形然后再次栅格化这些是非常奇怪的。我的印象是你让事情变得比需要的复杂得多。你也谈到圈子。哪些圈子？我猜你可能想围绕你的观点画圈，但这不是你在做什么。逐步弄清楚事情可能会有所帮助，首先弄清楚如何获得您想要的一般结果，然后如何使其以太平洋为中心。
以下是代码第一部分的清理版本。它还使它具有可重现性。您需要在代码中创建示例，如下所示：
lon <- c(-177.87567, -83.2136, -82.59392, -82.4106, -81.17555, -80.9577)
lat <- c(-24.715167, 14.4018, 9.589192, 9.49275, 7.19675, 8.8527)
radius <- c(10, 1, 1, 1, 5, 1)
site <- c("MIG", "OBS", "NES", "NES;BRE", "OBS", "NES")
sp_csv_data <- data.frame(lon, lat, radius, site)
## cleaned up projection definitions
rob_pacific <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
longlat <- "+proj=longlat +datum=WGS84"
##Converting to raster
# Creating a empty raster at 0.5° resolution
rs <- raster(res=0.5, crs=lonlat)
values(rs) <- 1:ncell(rs)
sp_raster <- rasterize(sp_csv_data[,1:2], rs, sp_csv_data[,3])
## makes no sense: sp_raster <- resample(sp_raster, rs, resample = "ngb")
#converting into an SpatialPolygonsDataframe
species_sp <- as(sp_raster, "SpatialPolygonsDataFrame")
## makes no sense: species_sp <- spTransform(sp_raster, CRS(longlat))