微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

将坐标、半径和站点类型数据从 .csv 转换为 R 中的栅格

如何解决将坐标、半径和站点类型数据从 .csv 转换为 R 中的栅格

我在将 .csv 文件转换为 R 中的栅格时遇到了一些问题...我的 .csv 文件包含坐标(长和纬度)、半径(以度为单位)和站点类型。我能够将坐标转换为栅格并能够使用 st_buffer() 绘制圆圈,但我面临两个问题:

  1. 我无法将圆圈转换为栅格...我尝试使用 rasterize()fasterize(),但两者都不起作用,我得到的是一个空栅格层
  2. 我似乎无法根据网站类型对坐标和圆圈进行分类

知道我可能做错了什么吗?以及如何对我的圈子进行分类? 提前致谢!

这是我使用的代码

> 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"

####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,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.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() 都给了我一个空的栅格层...

解决方法

rasterize 最终不工作的原因很明显:vector 和raster 的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,5,1)
site <- c("MIG","OBS","NES","NES;BRE","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,resample = "ngb")
    
#converting into an SpatialPolygonsDataframe
species_sp <- as(sp_raster,"SpatialPolygonsDataFrame")
## makes no sense: species_sp <- spTransform(sp_raster,CRS(longlat))

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。