In einarhjorleifsson/fishdown: Does not matter.

On stations

SMB

In 2017 the tow-index for tows having "tognumer" less than 40 have been fixed to strata. The code used in 2018 to generate the "official" index can be found under /u2/reikn/R/SurveyWork/SMB/stodvar.r

What follows is a tidy mimicry of that code

library(geo)
library(tidyverse)
library(mar)
con <- connect_mar()

Station strata

The Station strata are stored in ops$einarhj.SMB_INDEX_STRATA_AREA:

index.strata <- 
  tbl_mar(con, "ops$einarhj.SMB_INDEX_STRATA_AREA") %>% 
  # NOTE: should fix this nomenclature upstream
  select(index, stdnewstrata, stdoldstrata)
glimpse(index.strata)

Note that table contains also variable "area", standing for Bormicon area (in case one wants fixed Bormicon for each index).

Getting all stations and add strata (for all townumbers less than 40):

st <-
  les_stod(con) %>%
  left_join(les_syni(con)) %>% 
  filter(synaflokkur_nr == 30) %>% 
  #      NOTE: This should be fixed in fiskar.stodvar
  #      2019-02-13 Hoski to double check and then send in a correction to Sigrun
  mutate(tognumer = ifelse(synis_id == 474396, 1, tog_nr),
         index = reitur * 100 + tognumer) %>% 
  left_join(index.strata, by = "index")

Checking - all townumbers less than 40 should have a defined strata:

st %>% 
  mutate(less40 = ifelse(tognumer < 40, 1, 0),
         def.old = ifelse(!is.na(stdoldstrata), 1, 0),
         def.new = ifelse(!is.na(stdnewstrata), 1, 0)) %>% 
  group_by(less40) %>% 
  summarise(n = n(),
            def.old = sum(def.old, na.rm = TRUE),
            def.new = sum(def.new, na.rm = TRUE))

So all tows (19142) less than 40 have assigned a strata. And 123 tows are above 40 and they, as expected do not have a strata. To assign a strata to those station (not needed in any future calculations of indices) we use the geoinside approach and generate two new variables, oldstrata and newstrata (for all the stations). The values in these variables is then used in the stdoldstrata and stdnewstrata in cases townumber above 40.

# function in package husky changed so it returns a vector
inside_strata <- 
  function (stodvar, stratalist = ralllist, stratas = STRATAS) {
    stodvar$newstrata <- rep(0, nrow(stodvar))
    for (i in stratalist) {
      j <- geo::geoinside(stodvar, reg = stratas[[i]], robust = T, 
                          option = 0)
      if (length(j) > 0) 
        stodvar[j, "newstrata"] <- i
    }
    i <- stodvar$newstrata == 0
    if (any(i)) {
      i1 <- 1:length(i)
      i1 <- i1[i]
      cat("Strata fannst ekki fyrir stodvar", paste(i1, collapse = ","))
    }
    return(stodvar$newstrata)
  }


# Have to get the stations into R for the geoinside operation
st <-
  st %>% 
  collect(n = Inf) %>% 
  mutate(lon = kastad_lengd,
         lat = kastad_breidd) %>% 
  mutate(oldstrata = inside_strata(.,
                                   stratalist = husky::ralllist,
                                   stratas = husky::STRATAS),
         newstrata = inside_strata(.,
                                   stratalist = husky::smblist$nr,
                                   stratas = husky::NEWSTRATAS)) %>% 
  # Some fixing (should ideally not need this)
  mutate(oldstrata = ifelse(oldstrata %in% 23, 88, oldstrata),
         # townumber >= 40 are "aukastodvar" hence below is really not used in 
         # any official calculation  
         newstrata = ifelse(index %in% 71741 & ar %in% 2009, 37, newstrata),
         newstrata = ifelse(index %in% 61141 & ar %in% 2010, 42, newstrata)) %>% 
  mutate(oldstrata = ifelse(tognumer < 40, stdoldstrata, oldstrata),
         newstrata = ifelse(tognumer < 40, stdnewstrata, newstrata)) %>% 
  select(-c(stdoldstrata, stdnewstrata)) %>% 
  geo::inside.reg.bc(.)

Comparison with God

Strata allocation:

load("/u2/reikn/R/SurveyWork/SMB/Stations.rdata")
God <- 
  STODVAR.all %>% 
  select(synis_id = synis.id, index, oldstrata, newstrata) %>% 
  as_tibble() %>% 
  arrange(synis_id)
Tidy <- 
  st %>% 
  select(synis_id, index, oldstrata, newstrata) %>% 
  arrange(synis_id)
identical(God, Tidy)

Now check for the Bormicon area

God <-
  STODVAR.all %>% 
  filter(tognumer < 40) %>% 
  select(synis_id = synis.id, area) %>% 
  as_tibble() %>% 
  arrange(synis_id)
Tidy <- 
  st %>%
  filter(tognumer < 40) %>% 
  select(synis_id, area) %>% 
  arrange(synis_id)
identical(God, Tidy)

Note, we have index that falls in different Bormicon areas in different years:

st %>% 
  filter(tognumer < 40) %>% 
  select(index, area) %>% 
  distinct() %>% 
  group_by(index) %>% 
  count() %>% 
  filter(n > 1)

Question if we should not fix the Bormicon area for all index for tognumer < 4?? Is already availble in ops$einarhj.SMB_INDEX_STRATA_AREA

Bottom line

A function (without messing with townumber > 40)

get_smb_stations <- function(con) {

  lesa_stodvar(con) %>%
    filter(synaflokkur == 30) %>%
    #      NOTE: This should be fixed in fiskar.stodvar
    #      2019-02-13 Hoski to double check and then send in a correction to Sigrun
    mutate(tognumer = ifelse(synis_id == 474396, 1, tognumer),
           index = reitur * 100 + tognumer) %>% 
    left_join(tbl_mar(con, "ops$einarhj.SMB_INDEX_STRATA_AREA") %>% 
                rename(oldstrata = stdoldstrata,
                       newstrata = stdnewstrata), by = "index")

}