I am trying to summarize some statistics in the grid that I made, however something fails when I try to do it.
This is my data
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>, …
and I did this raster
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))
but I can not represent the data in the grid
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)
Please help me to find this solution!!
CodePudding user response:
So I just saw that you shifted the coordinates with st_shift_longitude()
and therefore your bounding box is:
Bounding box: xmin: 303.22 ymin: -61.43 xmax: 303.95 ymax: -60.78
Do you really need it? That doesn't match with your defined extent
e <- extent(-70,-40,-68,-60)
And a bbox for WGS84 is suppose to be at max c(-180, -90, 180, 90)
.
Also, on your plot you are not instructed ggplot2
to do anything with the values of catk
. Grid
and rc3
do not have anything from catk
, is the result
object.
Anyway, I give a try to your problem even though I don't have access to your dataset. I represent on each cell sum_cat
from result
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)
Created on 2022-03-23 by the reprex package (v2.0.1)