Отфильтруйте последние точки, используя шейп-файл в R - PullRequest
0 голосов
/ 19 сентября 2018

Я пытаюсь отфильтровать точки тропических циклонов с помощью шейп-файла:

  1. Ссылка на шейп-файл

Я делаюпроблема, аналогичная предыдущей:

Фильтр дорожек TC (точки широты) по указанному сеточному полю

Но в этом случае я использую неправильный шейп-файл (не поле).

Я хотел отфильтровать только точки с тем же идентификатором (столбец «SN» в приведенных ниже примерах данных), который прошел через шейп-файл.Таким образом, будут фильтроваться только уникальные TC.

Вот подмножество данных:

dat <- structure(list(SN = c(200608L, 200608L, 200608L, 200608L, 200608L, 
200608L, 200608L, 200608L, 200608L, 200608L, 200608L, 200608L, 
200610L, 200610L, 200610L, 200612L, 200612L, 200612L, 200612L, 
200612L, 200612L, 200612L, 200709L, 200709L, 200709L, 200709L, 
200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 
200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 
200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 200709L, 
200709L, 200709L, 200709L, 200709L, 200709L), CY = c(8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 10L, 10L, 10L, 12L, 12L, 
12L, 12L, 12L, 12L, 12L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L), Year = c(2006L, 2006L, 2006L, 2006L, 2006L, 
2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 
2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2007L, 
2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 
2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 
2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 2007L, 
2007L, 2007L), Month = c(8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L), Day = c(9L, 9L, 
9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 9L, 9L, 9L, 14L, 
14L, 14L, 15L, 15L, 15L, 15L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 
14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 
17L, 17L, 17L, 18L, 18L, 18L, 18L, 19L, 19L, 19L), Hour = c(0L, 
6L, 12L, 18L, 0L, 6L, 12L, 18L, 0L, 6L, 12L, 18L, 0L, 6L, 12L, 
6L, 12L, 18L, 0L, 6L, 12L, 18L, 6L, 12L, 18L, 0L, 6L, 12L, 18L, 
0L, 6L, 12L, 18L, 0L, 6L, 12L, 18L, 0L, 6L, 12L, 18L, 0L, 6L, 
12L, 18L, 0L, 6L, 12L, 18L, 0L, 6L, 12L), Lat = c(23.7, 24.7, 
25.3, 25.9, 26.4, 27, 27.2, 27.4, 27.7, 28.1, 28.5, 28.9, 22.8, 
22.9, 22.4, 18.7, 19.8, 20.2, 21, 22.4, 23.9, 26.9, 17.4, 17, 
16.7, 16.6, 16.5, 16.4, 16.3, 16.2, 16, 15.8, 15.5, 15.6, 15.9, 
16.1, 16.6, 17.2, 17.9, 18.7, 19.4, 20.2, 21, 21.9, 22.7, 23.4, 
24, 24.4, 24.8, 25.2, 25.6, 26), Lon = c(128.4, 126.9, 125.3, 
123.9, 122.5, 121.1, 119.9, 118.7, 117.4, 115.8, 115, 114.4, 
119.8, 118.6, 117.7, 131.3, 132.4, 133.9, 135.6, 137.2, 139.1, 
140.4, 135.4, 135.2, 134.8, 134.3, 133.7, 133, 132.3, 131.5, 
130.7, 129.9, 129.2, 128.5, 128, 127.6, 127.1, 126.6, 126.1, 
125.6, 124.8, 124.2, 123.5, 122.8, 122, 121.2, 119.9, 119.4, 
119, 118.6, 118.3, 118)), class = "data.frame", row.names = c(NA, 
-52L))

Полная ссылка данных здесь

ОжидаетсяРезультаты:

Мои ожидаемые результаты аналогичны ссылке выше.Два дополнительных столбца в файле csv, которые указывают, находятся ли TC в шейп-файле.Каждый TC имеет уникальный идентификатор: столбец SN.Таким образом, если точки-точки с одинаковым номером SN находятся в шейп-файле, они помечаются как ИСТИНА (все точки с одинаковым номером SN).

dat[SN %in% dat[inBounds == TRUE, unique(SN)], passesThroughBox := T ]
dat[is.na(passesThroughBox), passesThroughBox := F]

Любые предложения о том, как я могусделать это в R?

С уважением, Линдз

1 Ответ

0 голосов
/ 19 сентября 2018

Не лучшие выборочные данные, так как ни одно наблюдение, кажется, не проходит через bbox шейп-файла. Тем не менее, это мое решение проблемы.

library( tidyverse )
library( sf )

sf и решение tidyverse

#first, get the boundaries of the shapefile
bbox <- read_sf(dsn = "./Bicol_Region/Bicol_region.shp") %>% st_bbox()

#      xmin      ymin      xmax      ymax 
# 122.29929  11.71056 124.42607  14.50061 


dat %>% 
  #check if measurement is withing boundaries (yes, or not (no))
  mutate( passes_through_box = ifelse( Lat >= bbox[2] & Lat <= bbox[4] & Lon >= bbox[1] & Lon <= bbox[3], "yes", "no" ) ) %>%
  #group by SN
  group_by( SN ) %>%
  #check if any value of 'passes_through_box' in a group is "yes", if so 'yes', else 'no'
  mutate( passes_through_box_anyime = ifelse( any( passes_through_box == "yes"), "yes", "no" ) )

# # A tibble: 52 x 10
# # Groups:   SN [4]
#       SN    CY  Year Month   Day  Hour   Lat   Lon passes_through_box passes_through_box_anyime
#    <int> <int> <int> <int> <int> <int> <dbl> <dbl> <chr>              <chr>                    
# 1 200608     8  2006     8     9     0  23.7  128. no                 no                       
# 2 200608     8  2006     8     9     6  24.7  127. no                 no                       
# 3 200608     8  2006     8     9    12  25.3  125. no                 no                       
# 4 200608     8  2006     8     9    18  25.9  124. no                 no                       
# 5 200608     8  2006     8    10     0  26.4  122. no                 no                       
# 6 200608     8  2006     8    10     6  27    121. no                 no   

все sf решение

library( tidyverse )
library( sf )

#first, get the boundaries of the shapefile
bbox <- read_sf(dsn = "./Bicol_Region/Bicol_region.shp") %>% st_bbox() %>% st_as_sfc( crs = "+proj=longlat +datum=WGS84" )

#      xmin      ymin      xmax      ymax 
# 122.29929  11.71056 124.42607  14.50061 

#create aspatial data.frame
spdf <- st_as_sf( x = dat,
                  coords = c("Lon", "Lat"),
                  crs = "+proj=longlat +datum=WGS84" )

spdf %>% 
  #check if a line the spatial df intersecgt with the defined boundary-box
  mutate( passes_through_box = as.numeric( st_intersects(spdf, bbox) ) ) %>%
  group_by( SN ) %>%
  mutate( passes_through_box_anyime = ifelse( any( passes_through_box == 1), "yes", "no" ) )

для полноты картины, другой sf-подход

library( tidyverse )
library( sf )

#first, get the boundaries of the shapefile and create a box
bbox <- read_sf(dsn = "./Bicol_Region/Bicol_region.shp") %>% st_bbox() %>% st_as_sfc( crs = "+proj=longlat +datum=WGS84" )

#create a spatial points data.frame using the sample data provided
spdf <- st_as_sf( x = dat,
                  coords = c("Lon", "Lat"),
                  crs = "+proj=longlat +datum=WGS84" )

#create a spatial lines data.frame, bases on lat-lon groupes by SN
sldf <- spdf %>%
  group_by( SN ) %>%
  summarise( m = mean( Year ) ) %>%
  st_cast( "LINESTRING" ) %>%
  select( -m )

#let's see the printed results
library(mapview)
mapview( list(bbox, sldf) )

enter image description here

#any intersections? ... #nope, no intersections
st_intersects( bbox, sldf )

# although coordinates are longitude/latitude, st_intersects assumes that they are planar
# Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
# 1: (empty)
...