将坐标、半径和站点类型数据从 .csv 转换为 R 中的栅格 [英] Converting coordinates, radius and site type data from .csv to raster in R

问题描述

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

I am having a bit of a problem converting .csv files into raster in R... My .csv file contains coordinates (long and lat) radius (in deg) and site type. I was able to convert the coordinates into raster and was able to plot the circles using st_buffer() but I am facing two problems:

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

1. I can't convert the circles into a raster... I tried with rasterize() and fasterize() and both did not work all I'm getting is an empty raster layer
2. I can't seem to classify the coordinates and circles according to the site type

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

Any idea of what I might be doing wrong? and how can I classify my circles? Thank you in advance!

这是我使用的代码:

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

I want to convert sp_robinson2 into a raster called sp_robinson3 but as you can see both fasterize()and rasterize()are giving me an empty raster layer...

推荐答案

rasterize 最终不起作用的原因很明显:vector 和 raster 的 crs 不匹配.但是，您能否再编辑一下您的问题以解释您想要实现的目标?创建一个栅格，然后创建多边形，然后再次栅格化它们是非常奇怪的.我的印象是你让事情变得比需要的复杂得多.你也谈论圈子.哪些圈子?我猜你可能想要围绕你的观点圈起来，但这不是你在做什么.逐步弄清楚事情可能会有所帮助，首先弄清楚如何获得您想要的一般结果，然后如何使其以太平洋为中心.

The reason why rasterize does not work in the end is obvious: the crs of the vector and raster do not match. But can you edit your question a bit more to explain what you want to achieve? It is very odd to create a raster and then polygons and then rasterize these again. My impression is that you are making things much more complicated than need be. You also talk about circles. Which circles? I am guessing you may want circles around your points, but that is not what you are doing. It would probably be helpful to figure out things step by step, first figure out how to get the general result you want, then how to get it Pacific centered.

以下是代码第一部分的清理版本.这也使它具有可复制性.您需要在代码中创建示例，如下所示:

Below is a cleaned up version of the first part of your code. It also makes it reproducible. You need to create example in code, like this:

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 中的栅格的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！