如何在地图 sf 的网格中绘制我自己的数据，但 return 真空

Question

我试图在我制作的网格中总结一些统计数据，但是当我尝试这样做时出现了一些失败。

这是我的数据

head(catk)
Simple feature collection with 6 features and 40 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 303.22 ymin: -61.43 xmax: 303.95 ymax: -60.78
Geodetic CRS:  WGS 84
# A tibble: 6 × 41
     X1   day month  year c1_id greenweight_caught_kg obs_haul_id
  <dbl> <dbl> <dbl> <dbl> <dbl>                 <dbl> <lgl>      
1     1     4    12  1997 26529                  7260 NA         
2     2     4    12  1997 26530                  7920 NA         
3     3     4    12  1997 26531                  4692 NA         
4     4     4    12  1997 26532                  5896 NA         
5     5     4    12  1997 26533                    88 NA         
6     6     5    12  1997 26534                  7040 NA         
# … with 34 more variables: obs_logbook_id <lgl>, obs_haul_number <lgl>,
#   haul_number <dbl>, vessel_name <chr>, vessel_nationality_code <chr>,
#   fishing_purpose_code <chr>, season_ccamlr <dbl>,
#   target_species <chr>, asd_code <dbl>, trawl_technique <lgl>,
#   catchperiod_code <chr>, date_catchperiod_start <date>,
#   datetime_set_start <dttm>, datetime_set_end <dttm>,
#   datetime_haul_start <dttm>, datetime_haul_end <dttm>, …

我做了这个光栅

an <- getData("GADM", country = "ATA", level = 0)
an@data$NAME_0

e <- extent(-70,-40,-68,-60)
rc <- crop(an, e)
proj4string(rc) <- CRS("+init=epsg:4326")
rc3 <- st_as_sf(rc)

catk <- st_as_sf(catk, coords = c("Longitude", "Latitude"),  crs = 4326) %>%
  st_shift_longitude()

Grid <- rc3 %>%
  st_make_grid(cellsize = c(1,0.4)) %>% # para que quede cuadrada
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(cellid = row_number()) 

result <- Grid %>%
  st_join(catk) %>% 
  group_by(cellid) %>% 
  summarise(sum_cat = sum(Catcht))

但我无法表示网格中的数据

ggplot() +
  geom_sf(data = Grid, color="#d9d9d9", fill=NA) +
  geom_sf(data = rc3) +
  theme_bw() + 
  coord_sf() +
  scale_alpha(guide="none")+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides( colour = guide_legend()) +
  theme(panel.background = element_rect(fill = "#f7fbff"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  theme(legend.position = "right")+
  xlim(-69,-45)

fail plot

请帮我找到这个解决方案！！

Answer 1

所以我刚刚看到您用 st_shift_longitude() 移动了坐标，因此您的边界框是：

Bounding box:  xmin: 303.22 ymin: -61.43 xmax: 303.95 ymax: -60.78

你真的需要吗？这与您定义的范围不匹配

e <- extent(-70,-40,-68,-60)

WGS84 的 bbox 应该最大 c(-180, -90, 180, 90)。

此外，在您的绘图中，您没有被指示 ggplot2 使用 catk 的值做任何事情。 Grid 和 rc3 没有来自 catk 的任何东西，是 result 对象。

无论如何，即使我无法访问您的数据集，我也会尝试解决您的问题。我代表 result

中的每个单元格 sum_cat

library(raster)
library(sf)
library(dplyr)
library(ggplot2)

# Mock your data
catk <- structure(list(Longitude = c(-59.0860203764828, -50.1352159580643, 
                                     -53.7671292009259, -67.9105254106185, -67.5753491797446, -51.7045571975837, 
                                     -45.6203830411619, -61.2695183762776, -51.6287384188695, -52.244074640978, 
                                     -45.4625770258213, -51.0935832496694, -45.6375681312716, -44.744215508174, 
                                     -66.3625310507564), Latitude = c(-62.0038884948778, -65.307178606448, 
                                                                      -65.8980199769778, -60.4475595973147, -67.7543165093134, -60.4616894158005, 
                                                                      -67.9720957068844, -62.2184680275876, -66.2345680342004, -64.1523442367459, 
                                                                      -62.5435163857161, -65.9127866479611, -66.8874734854608, -62.0859917484373, 
                                                                      -66.8762861503705), Catcht = c(18L, 95L, 32L, 40L, 16L, 49L, 
                                                                                                     22L, 14L, 86L, 45L, 3L, 51L, 45L, 41L, 19L)), row.names = c(NA, 
                                                                                                                                                                 -15L), class = "data.frame")


# Start the analysis
an <- getData("GADM", country = "ATA", level = 0)
e <- extent(-70,-40,-68,-60)
rc <- crop(an, e)
proj4string(rc) <- CRS("+init=epsg:4326")

rc3 <- st_as_sf(rc) 

# Don't think you need st_shift_longitude, removed
catk <- st_as_sf(catk, coords = c("Longitude", "Latitude"),  crs = 4326) 


Grid <- rc3 %>%
  st_make_grid(cellsize = c(1,0.4)) %>% # para que quede cuadrada
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(cellid = row_number()) 




result <- Grid %>%
  st_join(catk) %>% 
  group_by(cellid) %>% 
  summarise(sum_cat = sum(Catcht))



ggplot() +
  geom_sf(data = Grid, color="#d9d9d9", fill=NA) +
  # Add here results with my mock data by grid
  geom_sf(data = result %>% filter(!is.na(sum_cat)), aes(fill=sum_cat)) +
  geom_sf(data = rc3) +
  theme_bw() + 
  coord_sf() +
  scale_alpha(guide="none")+
  xlab(expression(paste(Longitude^o,~'O'))) +
  ylab(expression(paste(Latitude^o,~'S')))+
  guides( colour = guide_legend()) +
  theme(panel.background = element_rect(fill = "#f7fbff"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  theme(legend.position = "right")+
  xlim(-69,-45)

^{由 reprex package (v2.0.1)}

于 2022-03-23 创建

如何在地图 sf 的网格中绘制我自己的数据，但 return 真空

How can plot my own data in a grid in a map sf but return vacum

gis

maps

shapes

ggplot2

sf