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

Time:03-25

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)

fail plot

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)

  • Related