GenerateOtherCatchesOnNodesTheHardLongWay.R

 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 ##!!find landings per month per species for each igraph node!!!##
 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 
 ## find landings per month per species for each igraph node
 ## to be used for the depletion of the pops each month that would result from
 ## the not explicitly described "other" activities
 ## in the simulation, e.g. from foreign countries...
 ## but should also be used for vessels not equipped with VMS i.e. smaller vessels...

 
 ## POP SPE
 ## IBM parametrisation
 ## Francois Bastardie (DTU-Aqua)
 ## outputs: mainly .dat files to be used for the IBM simulations

# GENERAL SETTINGS
  general <- list()
  general$main.path             <- file.path("C:", "Users", "fbas", "Documents", "GitHub", "DISPLACE_input_raw")
  general$main.path.code        <- file.path("C:", "Users", "fbas", "Documents", "GitHub", "DISPLACE_R_inputs")
  general$main_path_input       <- file.path("C:", "Users", "fbas", "Documents", "GitHub", "DISPLACE_input_myfish")

  general$igraph                <- 11
  general$case_study            <- "baltic_only"
  general$case_study_countries  <- c("DEN", "SWE", "DEU")    # for the Baltic only
  general$a.year                <- "2012"
  general$a.country             <- "DEN"
  #general$a.country             <- "DEU"
  #general$a.country             <- "SWE"


  general$igraph                <- 56
  general$case_study            <- "myfish"
  general$case_study_countries  <- c("DEN")    # for the Baltic only
  general$a.year                <- "2012"
  general$a.country             <- "DEN"


  general$region.of.interest   <- "BE_and_BW_and_Kattegat"
           
 
 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 ##!!!!!!!!!!!!!!!!!!!!!UTILS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##
 ##!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##

  library(maps)
 
  ICESarea2 <-
function (tacsat, string = TRUE)
{
    library(sp)
    ICES.area <- rep(NA, dim(tacsat)[1])
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(8.648336, 7.034822, 7.357525, 9.083985, 9.608377,
            10.22958, 10.689431, 11.084742, 11.617201, 12.068985,
            11.972174, 10.59262, 9.971417, 9.39862, 8.648336),
        pol.y = c(57.08073, 57.99182, 58.20964, 58.87187, 59.32015,
            59.86417, 59.99375, 59.8804, 58.96783, 58.0774, 57.4653,
            57.74247, 57.50441, 57.10708, 57.08073)) > 0] <- ifelse(string, 'kask', '0')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(10.59262, 11.97217, 12.15883, 12.70796, 13.12992,
            12.80622, 12.95073, 12.72185, 12.45127, 12.29556,
            12.13384, 11.99063, 11.58487, 11.58487, 11.63281,
            11.49492, 11.3094, 11.27652, 10.71374, 10.70218,
            10.24553, 10.19351, 10.42472, 10.59262), pol.y = c(57.74247,
            57.4653, 57.48032, 56.94085, 56.46389, 56.36135,
            56.19091, 56.16918, 56.29535, 56.12728, 55.49119,
            55.28764, 55.63113, 55.91101, 55.90623, 55.94866,
            55.97965, 56.00988, 56.14253, 56.25853, 56.49587,
            57.11107, 57.63566, 57.74247)) > 0] <-ifelse(string,'kask', '0')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(-4, 7, 8, 7, 7, 7.906163, -4, -4.6, -4.6, -4),
        pol.y = c(62, 62, 61.5, 60, 58, 57.5, 57.5, 57.3, 58.2,
            58.4)) > 0] <- ifelse(string,'nsea', '1')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(7.906163, 8.6488368, 8.5, 10.3, 10.3, 9.2,
            9.2, 7.11, -1, -4, -1.78), pol.y = c(57.5, 57.08073,
            57, 57.3, 57, 56.2, 52.7, 53.5, 53.5, 56.1, 57.5)) >
        0] <- ifelse(string,'nsea', '1')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(0, 0.5, 7.5, 7.5), pol.y = c(53.5, 51, 51,
            53.5)) > 0] <- ifelse(string,'nsea', '1')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(12, 12, 11.94, 11.97, 12, 12, 15, 15, 14.78,
            14.2, 13.7, 12.81, 12.44), pol.y = c(55.3, 54.75,
            54.67, 54.56, 54.56, 53.5, 53.5, 55, 55.3, 55.4,
            55.5, 55.38, 55.33)) > 0] <- ifelse(string,'2224', '2')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(14.2, 14.78, 15, 15, 18, 18, 14.2), pol.y = c(55.4,
            55.3, 55, 53, 53, 56.5, 56.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(18, 18, 22, 22), pol.y = c(56.5, 53.5, 53.5,
            56.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(18, 18, 18.32, 18.32, 19, 19, 16, 16), pol.y = c(56.5,
            57, 57, 57.5, 57.925, 59.762, 59.762, 56.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(19, 19, 18.45, 18.3, 18, 18, 21.5, 21.72, 21.98,
            22.17, 22.24, 21.93), pol.y = c(58.5, 57.9, 57.58,
            57, 57, 56.5, 56.5, 57.57, 57.97, 58.04, 58.15,
            58.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(21.5, 21.72, 21.98, 22.17, 22.24, 22.24, 23,
            25, 25), pol.y = c(56.5, 57.57, 57.97, 58.04, 58.15,
            58.35, 58.5, 58.5, 56.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(19, 17.975, 21.6, 21.8, 23.325, 23.325, 23.191,
            23, 23, 23.5, 23.6, 24, 23.692, 22.5, 22.1, 21.92,
            19), pol.y = c(59.762, 60.5, 60.5, 60.7, 60.5, 59.965,
            59.867, 59.827, 59, 59, 59.05, 58.75, 59.5, 59.5,
            58.35, 58.5, 58.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(16.5, 16.5, 19.7, 19.7, 22.6, 21.4), pol.y = c(60.5,
            63.7, 63.7, 63.5, 63.5, 60.5)) > 0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(19.7, 19.7, 25.7, 25.7, 19.7), pol.y = c(63.7,
            63.5, 63.5, 67, 67)) > 0] <-ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(23.325, 23.325, 23.191, 23, 23, 30.5, 30.5),
        pol.y = c(60.5, 59.965, 59.867, 59.827, 59, 59, 60.5)) >
        0] <- ifelse(string,'2532', '3')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(12.297, 12.13, 12.45, 12.81, 12.94, 13.21,
            12.5, 12.448), pol.y = c(56.13, 55.48, 55.31, 55.38,
            55.41, 55.71, 56.29, 56.305)) > 0] <- ifelse(string,'2224', '2')
    ICES.area[point.in.polygon(point.x = tacsat$SI_LONG, point.y = tacsat$SI_LATI,
        pol.x = c(10.1, 10.75, 10.71, 11.58, 11.58, 11.99, 11.94,
            11.97, 12, 12, 9.3, 9.3), pol.y = c(56.6, 56.3, 56.15,
            55.9, 55.65, 55, 54.67, 54.56, 54.56, 53.75, 53.75,
            56.6)) > 0] <- ifelse(string,'2224', '2')
    return(ICES.area)
}

 ICESrectangle <- function (dF) 
{
    rectChar1n2 <- as.integer(2 * (dF[, "SI_LATI"] - 35.5))
    rectChar3 <- ifelse(dF[, "SI_LONG"] <= -40, "A", ifelse(dF[, 
        "SI_LONG"] <= -30, "B", ifelse(dF[, "SI_LONG"] <= -20, 
        "C", ifelse(dF[, "SI_LONG"] <= -10, "D", ifelse(dF[, 
            "SI_LONG"] <= 0, "E", ifelse(dF[, "SI_LONG"] <= 10, 
            "F", ifelse(dF[, "SI_LONG"] <= 20, "G", ifelse(dF[, 
                "SI_LONG"] <= 30, "H", "I"))))))))
    rectChar4 <- as.integer(dF[, "SI_LONG"]%%10)
    rectID <- paste(rectChar1n2, rectChar3, rectChar4, sep = "")
    return(rectID)
}


  ICESrectangle2LonLat <- function (statsq, midpoint = F) {
    part1 <- substr(statsq, 1, 2)
    part2 <- substr(statsq, 3, 4)
    labels <- 0:90
    latlabels <- ifelse(labels < 10, paste("0", labels, sep = ""),
        as.character(labels))
    latvalues <- seq(35.5, 80.5, 0.5) + 0.25
    lonlabels <- paste(rep(LETTERS[2:8], rep(10, 7)), rep(0:9,
        7), sep = "")
    lonvalues <- (-40:29) + 0.5
    indx <- match(part1, latlabels)
    lat <- latvalues[indx]
    indx <- match(part2, lonlabels)
    lon <- lonvalues[indx]
    if (any(is.na(lat)) | any(is.na(lon)))
        warning("Some stat squares have not been recognised.")
    if (midpoint == F) {
        lat <- lat - 0.25
        lon <- lon - 0.5
    }
    return(data.frame(SI_LATI = lat, SI_LONG = lon))
   }


    c.listquote <- function (...)
   {
    args <- as.list(match.call()[-1])
    lstquote <- list(as.symbol("list"))
    for (i in args) {
        if (class(i) == "name" || (class(i) == "call" && i[[1]] !=
            "list")) {
            i <- eval(substitute(i), sys.frame(sys.parent()))
        }
        if (class(i) == "call" && i[[1]] == "list") {
            lstquote <- c(lstquote, as.list(i)[-1])
        }
        else if (class(i) == "character") {
            for (chr in i) {
                lstquote <- c(lstquote, list(parse(text = chr)[[1]]))
            }
        }
        else stop(paste("[", deparse(substitute(i)), "] Unknown class [",
            class(i), "] or is not a list()", sep = ""))
    }
    return(as.call(lstquote))
   }

 # function to assign a ICES area code to each ping
# based on a simplified version of the ICES_areas shape file
# compatible area code when uploading to fishframe


 longlat_to_ICESareas <- function(obj){

    library(sp)

    ICES.area <- rep(NA, dim(obj)[1])  # init

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(68.5,30.7,26,26,30,30,68.5),
              pol.y=c(63,63,70.648,72,72,90,90))>0] <- "Ib"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(33.7,34.55,35.28,36.38,37.57,
              38.31,39.05,39.61,41.24,42.81,43.06,44.68,43.51,43.18,41.73,41.56,40.66,40.51,39.76,38.96,37.74,36.61,35.70,33.70),
              pol.y=c(73.98,74.18,74.36,74.71,75.14,75.45,75.84,76.26,76.61,76.90,76.9,76.75,75.99,75.39,74.82,73.98,73.17,72.20,
              72.26,72.62,73.04,73.37,73.56,73.98))>0] <- "Ia"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-0.2,7.21,7.28,7.83,8.65,9.33,9.83,
              10.29,9.94,9.7,8.75,7.93,7.42,6.73,5.64,5.01,4.74,4.32,4,3.73,3.57,3.4,3.27,3.19,3.16,3.15,3.18,3.24,3.31,
              3.42,3.27,3.18,2.79,2.24,1.79,1.44,1.26,0.72,0.04,-0.489,-1.309,-1.559,-2.169,-2.539,-3.189,-3.729,-4.189,
              -4.559,-5.579,-5.599,-5.669,-5.779,-6.249,-6.619,-5.329,-4.189,-3.419,-2.389,-1.559,-0.609,0.08,0.68,1.18,
              1.46,1.72,1.94,2.09,2.25,2.35,2.39,2.38,2.31,2.22,2.06,1.89,1.68,1.48,1.08,0.34,-0.199), pol.y=c(73.5,73.5,
              73.45,73.14,72.76,72.49,72.31,72.18,71.98,71.91,71.64,71.36,71.13,70.79,70.17,69.79,69.56,69.32,69.1,68.86,
              68.69,68.46,68.23,67.98,67.77,67.57,67.37,67.18,67.01,66.84,66.43,66.39,66.23,65.95,65.64,65.38,65.32,65.08,
              64.72,64.43,64.84,64.92,65.13,65.22,65.39,65.47,65.55,65.59,65.69,65.96,66.22,66.47,67.09,67.61,67.77,67.96,
              68.1,68.33,68.55,68.86,69.14,69.44,69.76,69.97,70.21,70.43,70.63,70.89,71.14,71.35,71.61,71.83,72.01,72.24,
              72.43,72.6,72.75,72.99,73.31,73.5))>0] <- "IIa1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-11,-0.2,0.34,1.08,1.48,1.68,
              1.89,2.06,2.22,2.31,2.38,2.39,2.35,2.25,2.09,1.94,1.72,1.46,1.18,0.68,0.08,-0.609,-1.559,-2.389,-3.419,-4.189,
              -5.329,-6.619,-6.249,-5.779,-5.669,-5.599,-5.579,-4.559,-4.189,-3.729,-3.189,-2.539,-2.169,-1.559,-1.309,
              -0.489,0.04,0.72,1.26,1.44,1.79,2.24,2.79,3.18,3.27,3.42,3.31,3.24,3.18,3.15,3.16,3.19,3.27,3.4,3.57,3.73,
              4,4.32,4.74,5.01,5.64,6.73,7.42,7.93,8.75,9.7,9.94,10.29,9.83,9.33,8.65,7.83,7.28,7.21,30,30,26,26,11,-4,-4),
              pol.y=c(63,73.5,73.5,73.31,72.99,72.75,72.6,72.43,72.24,72.01,71.83,71.61,71.35,71.14,70.89,70.63,70.43,70.21,
              69.97,69.76,69.44,69.14,68.86,68.55,68.33,68.1,67.96,67.77,67.61,67.09,66.47,66.22,65.96,65.69,65.59,65.55,
              65.47,65.39,65.22,65.13,64.92,64.84,64.43,64.72,65.08,65.32,65.38,65.64,65.95,66.23,66.39,66.43,66.84,67.01,
              67.18,67.37,67.57,67.77,67.98,68.23,68.46,68.69,68.86,69.1,69.32,69.56,69.79,70.17,70.79,71.13,71.36,71.64,
              71.91,71.98,72.18,72.31,72.49,72.76,73.14,73.45,73.5,73.5,72,72,69.8,62,62,63))>0] <- "IIa2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(7.21,-0.2,-0.48,-1.88,-2.7,-5,-4.38,
             -4.29,-4.19,-4.3,-4.09,-2.52,-2.1,-1.6,0.8,1.12,1.71,3.06,4.07,4.55,5.19,6.39,6.51,6.74,7.06,7.21), pol.y=c(73.5,
             73.5,73.6,73.94,74.09,74.21,74.5,75,75.3,76.05,76.18,76.57,76.67,76.56,76,75.87,75.64,75.21,74.96,74.86,74.69,74.34,
             74.13,73.89,73.6,73.5))>0]<- "IIb1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-0.199,-0.479,-1.879,-2.699,-4.999,
             -4.379,-4.289,-4.189,-4.299,-4.089,-2.519,-2.099,-1.599,0.8,1.12,1.71,3.06,4.07,4.55,5.19,6.39,6.51,6.74,7.06,7.21,
             30,30,-11), pol.y=c(73.5,73.5,73.6,73.94,74.09,74.21,74.5,75,75.3,76.05,76.18,76.57,76.67,76.56,76,75.87,75.64,75.21,
             74.96,74.86,74.69,74.34,74.13,73.89,73.6,73.5,73.5,90,90))>0] <- "IIb2"

    #ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(7,7,8,8,8.5,10.3,10.3,10.1,10.75,10.71,11.58,
    #          11.58,12,12.31,12.47,13,13,9), pol.y=c(58.3,57.5,57.5,57,57,57.3,57,56.6,56.3,56.15,55.9,55.65,55.56,56.13,56.3,56,
    #          60.5,60.5))>0] <- "IIIa"

     ICES.area[point.in.polygon(point.x = obj[,"x"], point.y = obj[,"y"], 
        pol.x = c(8.648336, 7.034822, 7.357525, 9.083985, 9.608377, 
            10.22958, 10.689431, 11.084742, 11.617201, 12.068985, 
            11.972174, 10.59262, 9.971417, 9.39862, 8.648336), 
        pol.y = c(57.08073, 57.99182, 58.20964, 58.87187, 59.32015, 
            59.86417, 59.99375, 59.8804, 58.96783, 58.0774, 57.4653, 
            57.74247, 57.50441, 57.10708, 57.08073)) > 0] <-  "IIIan"
                                 
     ICES.area[point.in.polygon(point.x = obj[,"x"], point.y = obj[,"y"], 
        pol.x = c(10.59262, 11.97217, 12.15883, 12.70796, 13.12992, 
            12.80622, 12.95073, 12.72185, 12.45127, 12.29556, 
            12.13384, 11.99063, 11.58487, 11.58487, 11.63281, 
            11.49492, 11.3094, 11.27652, 10.71374, 10.70218, 
            10.24553, 10.19351, 10.42472, 10.59262), pol.y = c(57.74247, 
            57.4653, 57.48032, 56.94085, 56.46389, 56.36135, 
            56.19091, 56.16918, 56.29535, 56.12728, 55.49119, 
            55.28764, 55.63113, 55.91101, 55.90623, 55.94866, 
            55.97965, 56.00988, 56.14253, 56.25853, 56.49587, 
            57.11107, 57.63566, 57.74247)) > 0] <- "IIIas"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-4,7,8,7,7,-4,-4.6,-4.6,-4), pol.y=c(62,62,
              61.5,60,57.5,57.5,57.3,58.2,58.4))>0] <- "IVa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(8,8,8.5,10.3,10.3,9.2,9.2,7.11,-1,-4,-1.78),
              pol.y=c(57.5,57,57,57.3,57,56.2,52.7,53.5,53.5,56.1,57.5))>0] <- "IVb"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(0,0.5,7.5,7.5),
                pol.y=c(53.5,51,51,53.5))>0] <- "IVc"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-24,-27,-27,-24,-24),
                pol.y=c(62,62,63,63,62)) >0] <- "Va1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-27,-11,-11,-15,-15,-24,-24,-27),
                pol.y=c(68,68,63,63,62,62,63,63))>0] <- "Va2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-13.5,-15,-15,-14,-13.29,-13.5),
                pol.y=c(60,60,60,60.7,60.15,60)) > 0] <- "Vb1a"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-4,-5,-5,-8,-8,-7.5,-7.5,-8,-10,-10,-12,
                -13.5,-13.29,-13.99,-15,-15,-15,-11,-4), pol.y=c(60.5,60.5,60,60,60.5,60.5,61.25,61.5,61.5,60,60,60,60.15,60.71,
                60.49,62,63,63,63)) > 0] <- "Vb1b"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-8,-10,-10,-8,-7.5,-7.5,-8,-8), pol.y=c(60,
                60,61.5,61.5,61.25,60.5,60.5,60)) > 0] <- "Vb2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-4,-4,-4.6,-4.6,-4,-6,-8,-12,-12,-5,-5),
                pol.y=c(60.5,58.4,58.2,57.3,55,55,54.5,54.5,60,60,60.5)) > 0] <- "VIa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-15.99,-18,-18,-15,-13.5,-13.99,-14.57,
               -14.79,-14.88,-14.63,-14.34,-14.44,-14.54,-14.62,-14.72,-14.8,-14.89,-14.97,-15.04,-15.11,-15.19,-15.27,-15.34,-15.41,-15.47,
                -15.54,-15.6,-15.65,-15.7,-15.75,-15.79,-15.83,-15.87,-15.9,-15.92,-15.95,-15.97,-15.99,-15.99), pol.y = c(54.5,54.5,60,60,
                60,59.65,59.01,58.51,57.87,57.01,56.57,56.5,56.44,56.37,56.31,56.24,56.17,56.09,56.02,55.95,55.88,55.8,55.73,55.65,55.57,
                55.5,55.42,55.34,55.26,55.18,55.09,55.01,54.93,54.84,54.76,54.68,54.59,54.51,54.5)) >0] <- "VIb1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-12,-15.99,-15.99,-15.97,-15.95,-15.92,
                -15.9,-15.87,-15.83,-15.79,-15.75,-15.7,-15.65,-15.6,-15.54,-15.47,-15.41,-15.34,-15.27,-15.19,-15.11,-15.04,-14.97,-14.89,
                -14.8,-14.72,-14.62,-14.54,-14.44,-14.34,-14.63,-14.88,-14.79,-14.57,-13.99,-13.5,-12,-12), pol.y=c(54.5,54.5,54.51,54.59,
                54.68,54.76,54.84,54.93,55.01,55.09,55.18,55.26,55.34,55.42,55.5,55.57,55.65,55.73,55.8,55.88,55.95,56.02,56.09,56.17,56.24,
                56.31,56.37,56.44,56.5,56.57,57.01,57.87,58.51,59.01,59.65,60,60,54.5)) > 0] <- "VIb2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-6.3,-7.8,-2.5,-2.5), pol.y=c(55,52,52,55)) > 0] <- "VIIa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-8,-8,-12,-12),pol.y=c(54.5,52.5,52.5,54.5)) >0] <- "VIIb"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-16.06,-18,-18,-15.99,-15.99,-16,-16.01,-16.01,
                -16.01,-16,-15.99,-15.97,-15.96,-15.94,-15.91,-15.9,-15.89,-15.88,-15.86,-15.84,-15.88,-15.92,-15.95,-15.98,-16,-16.02,
                -16.04,-16.06,-16.06), pol.y=c(52.5,52.5,54.5,54.5,54.42,54.34,54.25,54.17,54.08,53.99,53.91,53.82,53.74,53.66,53.57,
                53.49,53.42,53.34,53.2,53.18,53.10,53.02,52.94,52.86,52.77,52.69,52.61,52.52,52.5)) >0] <- "VIIc1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-12,-16.06,-16.06,-16.04,-16.02,-16,-15.98,
                -15.95,-15.92,-15.88,-15.84,-15.86,-15.88,-15.89,-15.9,-15.91,-15.94,-15.96,-15.97,-15.99,-16,-16.01,-16.01,-16.01,
                -16,-15.99,-15.99,-12,-12), pol.y=c(52.5,52.5,52.52,52.61,52.69,52.77,52.86,52.94,53.02,53.1,53.18,53.26,53.34,53.42,
                53.49,53.57,53.66,53.74,53.82,53.91,53.99,54.08,54.17,54.25,54.34,54.42,54.5,54.5,52.5)) >0]<- "VIIc2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(2,-2,-2,-1.956,-1,2), pol.y=c(51,51,50.6,
                49.705,49,49)) >0] <- "VIId"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-2,-2.22,-5.17,-5.24,-7,-7,-5,-5,-1,-1,-1.956),
                pol.y=c(50.6,50.88,50.21,50,50,49.5,49.5,48,48,49.3,49.705)) > 0] <- "VIIe"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-7,-7,-5.2,-5.2,-2.4,-2.4,-5,-5,-6,-6),
                pol.y=c(50.5,50,50,50.25,51.4,51.9,51.9,51,51,50.5)) >0] <- "VIIf"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-9,-7,-7,-6,-6,-5,-5,-9),
                pol.y=c(50,50,50.5,50.5,51,51,52,52)) > 0] <- "VIIg"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-5,-8,-9,-9,-7,-7,-5,-5), pol.y=c(48,48,
                48,50,50,49,49,48)) >0] <- "VIIh"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-10.64,-11,-12,-12,-11.99,-11.87,-11.75,-11.64,
                -11.52,-11.39,-11.27,-11.14,-11.02,-10.89,-10.77,-10.68,-10.64), pol.y=c(48,48,48,48.43,48.42,48.39,48.36,48.33,48.3,
                48.27,48.25,48.23,48.21,48.19,48.17,48.03,48)) >0] <- "VIIj1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-9,-9,-10.64,-10.68,-10.77,-10.89,-11.02,-11.14,
                -11.27,-11.39,-11.52,-11.64,-11.75,-11.87,-11.99,-12,-12), pol.y=c(52.5,48,48,48.03,48.17,48.19,48.21,48.23,48.25,48.27,
                48.3,48.33,48.36,48.39,48.42,48.43,52.5)) >0] <- "VIIj2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-12,-12,-18,-18,-16.06,-16.07,-16.08,-16.09,-16.09,
                -16.09,-16.08,-16.07,-16.07,-16.07,-16.04,-16.02,-16,-15.96,-15.93,-15.9,-15.86,-15.82,-15.77,-15.73,-15.68,-15.63,-15.57,
                -15.52,-15.47,-15.42,-15.36,-15.3,-15.24,-15.17,-15.11,-15.04,-14.97,-14.89,-14.82,-14.74,-14.65,-14.57,-14.48,-14.39,-14.3,
                -14.22,-14.13,-14.04,-13.95,-13.86,-13.77,-13.67,-13.57,-13.47,-13.37,-13.27,-13.17,-13.07,-12.96,-12.85,-12.74,-12.64,-12.54,
                -12.43,-12.32,-12.22,-12.11,-12), pol.y=c(48.43,48,48,52.5,52.5,52.44,52.36,52.27,52.19,52.11,52.02,51.94,51.85,51.77,51.68,
                51.6,51.52,51.43,51.34,51.27,51.18,51.1,51.02,50.94,50.86,50.78,50.7,50.62,50.54,50.47,50.39,50.32,50.24,50.17,50.1,50.03,
                49.96,49.89,49.82,49.75,49.69,49.62,49.56,49.5,49.44,49.38,49.32,49.27,49.21,49.15,49.1,49.05,49,48.95,48.9,48.86,48.81,
                48.77,48.73,48.69,48.65,48.62,48.58,48.55,48.52,48.49,48.46,48.43)) >0] <- "VIIk1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-12,-12.11,-12.22,-12.32,-12.43,-12.54,-12.64,-12.74,
                -12.85,-12.96,-13.07,-13.17,-13.27,-13.37,-13.47,-13.57,-13.67,-13.77,-13.86,-13.95,-14.04,-14.13,-14.22,-14.3,-14.39,-14.48,
                -14.57,-14.65,-14.74,-14.82,-14.89,-14.97,-15.04,-15.11,-15.17,-15.24,-15.3,-15.36,-15.42,-15.47,-15.52,-15.57,-15.63,-15.68,
                -15.73,-15.77,-15.82,-15.86,-15.9,-15.93,-15.96,-15.99,-16.02,-16.04,-16.05,-16.07,-16.07,-16.08,-16.09,-16.09,-16.09,-16.08,
                -16.07,-16.06,-12,-12), pol.y=c(48.43,48.46,48.49,48.52,48.55,48.58,48.62,48.65,48.69,48.73,48.77,48.81,48.86,48.9,48.95,49,
                49.05,49.1,49.15,49.21,49.27,49.32,49.38,49.44,49.5,49.56,49.62,49.69,49.75,49.82,49.89,49.96,50.03,50.1,50.17,50.24,50.32,
                50.39,50.47,50.54,50.62,50.7,50.78,50.86,50.94,51.02,51.1,51.18,51.27,51.34,51.43,51.52,51.6,51.68,51.77,51.85,51.94,52.02,
                52.11,52.19,52.27,52.36,52.44,52.5,52.5,48.43)) >0] <- "VIIk2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-8,-8,-6,-6,-5,-5,-0.5,-0.5),
                pol.y=c(48,47.5,47.5,47,57,46,46,48)) >0] <- "VIIIa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-4,-4,-3,-3,-2,-2,0,0), pol.y=c(46,45.5,45.5,
                44.5,44.5,43,43,46)) >0] <- "VIIIb"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-11,-2,-2), pol.y=c(43,44.5,44.5,43)) >0] <- "VIIc"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-11,-10.64,-10.37,-9.89,-9.62,-10.95,-11),
                pol.y = c(46.32,48,48,47.77,47.45,46.88,46.34,46.32)) > 0] <- "VIIId1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-4,-4,-3,-3,-11,-11,-10.95,-9.62,-9.89,-10.37,-10.64,
                -9,-8,-8,-6,-5,-4,-4,-4), pol.y=c(46,45.5,45.5,44.5,44.5,46.32,46.34,46.88,47.45,47.77,48,48,48,47.5,47.5,47,47,46,46))>0] <- "VIIId2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(11,-13.31,-13.49,-13.8,-18,-18,-12,-11,-11),
                pol.y=c(46.32,44.72,44.07,43,43,48,48,48,46.32)) >0] <- "VIIIe1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-11,-13.8,-13.49,-13.31,-11,-11), pol.y=c(
                44.5,43,43,44.07,44.72,46.32,44.5)) >0] <- "VIIIe2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-11,-5.5,-5.5), pol.y=c(36,43,43,36))>0] <- "IXa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-17,-18,-13.8,-13.84,-13.64,-13.27,-13.27,-13.49,
                -13.78,-13.69,-12.73,-15.3,-17.9,-18), pol.y=c(36,43,43,42.88,42.04,41.38,41.13,40.06,38.75,38.17,36.03,36.04,36.02,36))
                >0] <- "IXb1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-18,-17.9,-15.3,-12.73,-13.69,-13.78,-13.49,
                -13.27,-13.27,-13.64,-13.84,-13.8,-11,-11), pol.y=c(36,36,36,36.04,36.03,38.17,38.75,40.06,41.13,41.38,42.04,42.88,43,43,
                36))>0] <- "IXb2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-18,-22.25,-20.62,-21.32,-23.91,-24.65,-25.79,
                -28.45,-29.95,-35.11,-35.26,-35.48,-31.76,-32.03,-42,-42,-18,-18), pol.y=c(36,36,37.58,39.16,40.97,41.35,41.91,42.34,
                42.05,41.02,40.04,38.74,36.03,36,36,43,43,36))>0] <- "Xa1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-22.25,-31.76,-32.03,-35.48,-35.26,-35.11,-29.95,
    -28.45,-25.79,-24.65,-23.91,-21.32,-20.62,-22.25), pol.y=c(36,36,36.03,38.74,40.04,41.02,42.05,42.34,41.91,41.35,40.97,39.16,37.58,
    36))>0] <- "Xa2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-18,-42,-42,-18,-18), pol.y=c(43,43,48,48,43))>0] <- "Xb"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-20.55,-24,-24,-18,-18,-42,-42,-41.5,-41,-40.5,
                -40,-39.5,-39,-38.5,-38.25,-38,-37.5,-37.2,-37,-36.77,-27,-27,-26.46,-25.09,-23.96,-23.27,-21.77,-20.57), pol.y=c(60,60,54.5,
                54.5,52.5,52.5,56.55,56.64,56.75,56.88,57.03,57.2,57.37,57.62,57.78,57.97,58.26,58.5,58.63,59,59,60.85,60.69,60.45,60.37,60.22,
                60.02,60))>0] <- "XIIa1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-20.55,-18.65,-17.32,-15.22,-15,-15,-18,-20.55),
                pol.y=c(60,60.05,60.11,60.44,60.49,60,60,60))>0] <- "XIIa2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-42,-36.77,-37,-37.2,-37.5,-38,-38.25,-38.5,-39,
                -39.5,-40,-40.5,-41,-41.5,-42,-42), pol.y=c(59,59,58.63,58.5,58.26,57.97,57.78,57.62,57.37,57.2,57.03,56.88,56.75,56.64,
                56.55,59))>0]<-"XIIa3"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-27,-27,-24,-15,-15,-15.22,-17.32,-18.65,-20.55,
                -21.76,-23.27,-23.96,-25.09,-26.46,-27), pol.y=c(60.85,62,62,62,60.49,60.44,60.11,60.05,60,60.02,60.22,60.37,60.45,
                60.69,60.85))>0]<-"XIIa4"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-18,-24,-24,-20.55,-18,-18), pol.y=c(54.5,54.5,
                60,60,60,54.5))>0]<-"XIIb"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-18,-42,-42,-18,-18), pol.y=c(48,48,52.5,52.5,
                48))>0]<-"XIIc"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-11,-27,-27,-40,-40,-11), pol.y=c(68,68,68.57,
                68.57,90,90))>0]<-"XIVa"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-27,-27,-36.77,-36.5,-36.35,-36.16,-35.96,-35.76,
                -35.5,-35.37,-35.15,-34.97,-34.65,-34.5,-34.31,-34,-33.7,-33.53,-33.27,-33,-32.5,-32.3,-32,-31.5,-31,-30.86,-30.61,-29.87,
                -29.25,-28.61,-27.69,-27), pol.y=c(60.85,59,59,59.35,59.5,59.75,60,60.25,60.55,60.75,61,61.25,61.5,61.6,61.75,61.98,62.25,
                62.45,62.5,62.56,62.69,62.75,62.87,63.03,63.25,63.31,63,62.23,61.79,61.44,61.06,60.85))>0]<-"XIVb1"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(-27,-27,-27.69,-28.61,-29.25,-29.87,-30.61,-30.86,
                -31,-31.5,-32,-32.3,-32.5,-33,-33.27,-33.53,-33.7,-34,-34.31,-34.5,-34.65,-34.97,-35.15,-35.37,-35.5,-35.76,-35.96,-36.16,
                -36.35,-36.5,-36.77,-42,-44,-44), pol.y=c(68.6,60.85,61.06,61.44,61.79,62.23,63,63.31,63.25,63.03,62.87,62.75,62.69,62.56,
                62.5,62.45,62.25,61.98,61.75,61.6,61.5,61.25,61,60.75,60.55,60.25,60,59.75,59.5,59.35,59,59,59,68.6))>0]<-"XIVb2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(12,12,11.94,11.97,12,12,15,15,14.78,14.2,13.7,
                12.81,12.44), pol.y=c(55.3,54.75,54.67,54.56,54.56,53.5,53.5,55,55.3,55.4,55.5,55.38,55.33))>0]<-"24"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(14.2,14.78,15,15,18,18,14.2), pol.y=c(55.4,55.3,
                55,53.7,53.7,56.5,56.5))>0]<-"25"  # corrected fba 14/12/2011

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(18,18,22,22),pol.y=c(56.5,53.5,53.5,56.5))>0]<-"26"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(18, 18, 18.32, 18.32, 19, 19, 16, 16), pol.y=c(56.5,
                57, 57, 57.5, 57.925, 59.762, 59.762, 56.5))>0]<-"27" # corrected fba 14/12/2011
                                                      
            
    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x = c(19, 19, 18.45, 18.3, 18, 18, 21.5, 21.72, 21.98,
            22.17, 22.24, 21.93), pol.y = c(58.5, 57.9, 57.58, 57, 57, 56.5, 56.5, 57.57, 57.97, 58.04, 58.15, 58.5))>0]<-"28-1"  # corrected fba 14/12/2011

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(21.5,21.72,21.98,22.17,22.24,22.24,23,25,25),
                pol.y=c(56.5,57.57,57.97,58.04,58.15,58.35,58.5,58.5,56.5))>0] <- "28-2"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(19,17.975,21.6,21.8,23.325,23.325,23.191,23,
                23,23.5,23.6,24,23.692,22.5,22.1,21.92,19), pol.y=c(59.762,60.5,60.5,60.7,60.5,59.965,59.867,59.827,59,59,59.05,
                58.75,59.5,59.5,58.35,58.5,58.5))>0]<-"29"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(16.5,16.5,19.7,19.7,22.6,21.4), pol.y=c(60.5,
                63.7,63.7,63.5,63.5,60.5))>0]<-"30"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(19.7,19.7,25.7,25.7,19.7), pol.y=c(63.7,63.5,
                63.5,67,67))>0]<-"31"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(23.325,23.325,23.191,23,23,30.5,30.5),
                pol.y=c(60.5,59.965,59.867,59.827,59,59,60.5))>0]<-"32"


    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(12.297,12.13,12.45,12.81,12.94,13.21,12.5,12.448),
                pol.y=c(56.13,55.48,55.31,55.38,55.41,55.71,56.29,56.305))>0]<-"23"

    ICES.area[point.in.polygon(point.x=obj[,"x"], point.y=obj[,"y"], pol.x=c(10.1,10.75,10.71,11.58,11.58,11.99,11.94,11.97,
                12,12,9.3,9.3), pol.y=c(56.6,56.3,56.15,55.9,55.65,55,54.67,54.56,54.56,53.75,53.75,56.6))>0]<-"22"



 return(ICES.area)
 }


 ##---------------------------------------
 ## READ IGRAPH CODE SQUARE---------------
 ##---------------------------------------
  #code_square_per_node <- read.table(file=file.path(general$main.path,  "igraph", 
  #                            paste("code_square_for_graph",general$igraph,"_points.dat", sep='')), header=TRUE) # built in ibm_param_vessels.r

  # replaced by:
  coord <- read.table(file=file.path(general$main_path_input, "graphsspe", paste("coord", general$igraph, ".dat", sep=""))) # build from the c++ gui
  coord <- as.matrix(as.vector(coord))
  coord <- matrix(coord, ncol=3)
  colnames(coord) <- c('SI_LONG', 'SI_LATI', 'idx.port')
  code_square_per_node <- cbind.data.frame(coord, code_square= ICESrectangle (coord))
  #plot(coord[,1], coord[,2])

  
  code_square_per_node$pt_graph <- 0:(nrow(code_square_per_node)-1)




  ##---------------------------------------
  ## IMPORT LANDINGS IN KILO DATA----------
  ##---------------------------------------
   cat(paste("import landings data from STECF (annual landings in tons per rectangle)\n",sep=''))

   # http://stecf.jrc.ec.europa.eu/data-reports
   # All countries from 2012 STECF effort regime working group---------     
   # annual landings only...NOT BY QUARTER (OR AT LEAST SEMESTER)!!
   landings_all <- read.table(file=file.path(general$main.path, "STECF",
                           'stecf_landings_per_rectangle_per_country_2012.csv'), sep=",", header=TRUE)

   # keep only the essentials...
   landings_all <-  landings_all[,c('COUNTRY','REG.GEAR.COD', 'VESSEL_LENGTH', 'RECTANGLE', 'SPECIES', 'X2012') ]
                            ##=> in TONS!
  
   #...and aggregate
   nm         <- names(landings_all)
   idx.col    <- grep('X2012', nm) # index columns with species weight
   library(data.table)
   DT  <- data.table(landings_all) # library data.table for fast grouping replacing aggregate()
   # AGGREGATE WEIGHT PER SPECIES
   eq1                <- c.listquote( paste ("sum(",nm[idx.col],",na.rm=TRUE)",sep="") )
   DT$X2010           <- as.double(DT$X2012)
   DT$REG.GEAR.COD    <- as.factor(DT$REG.GEAR.COD)
   DT$VESSEL_LENGTH   <- as.factor(DT$VESSEL_LENGTH)
   DT$RECTANGLE       <- as.factor(DT$RECTANGLE)
   DT$SPECIES         <- as.factor(DT$SPECIES)
   DT$COUNTRY         <- as.factor(DT$COUNTRY)
   land.agg           <- DT[,eval(eq1), by=list(COUNTRY, REG.GEAR.COD, VESSEL_LENGTH, RECTANGLE, SPECIES)]
   land.agg           <- data.frame(land.agg)
   colnames(land.agg) <- c("country","metier", "vsize", "rectangle", "species", "landings")
  
   # correct some levels....baltic_only_2010
   #levels(land.agg$vsize)
   #"NONE"    "O10M"    "o10t12m" "O10T12M" "o10t15m" "O10T15M" "o12t18m" "O12T18M" "o15m"    "O15M"    "o18t24m" "O18T24M" "o24t40m" "O24T40M" "o40m"   
   #"O40M"    "o8t10m"  "O8T10M"  "u10m"    "U10M"    "u8m"     "U8M" 
   #levels(land.agg$vsize) <- c("u12m", "u12m", "o10t12m","o10t12m","o10t15m","o10t15m","o12t18m","o12t18m","o15m", "o15m","o18t24m","o18t24m","o24t40m","o24t40m","o40m",   
   #                                   "o40m","o8t10m","o8t10m" ,"u10m","u10m","u8m","u8m") 
   
   # correct some levels....
   levels(land.agg$vsize)
   # "none"    "NONE"    "o10m"    "O10M"    "o10t12m" "O10T12M" "o10t15m" "O10T15M" "o12t18m" "O12T18M" "o15m"    "O15M"    "o18t24m" "O18T24M" "o24t40m" "O24T40M" "o40m"    "O40M"    "o8t10m"  "O8T10M"  "u10m"    "U10M"    "u8m"    
   # "U8M"
   levels(land.agg$vsize) <- c("u12m", "u12m", "o10t12m", "o10t12m", "o10t12m", "o10t12m", "o10t15m", "o10t15m", "o12t18m", 
                                "o12t18m", "o15m", "o15m", "o18t24m", "o18t24m", "o24t40m", "o24t40m", "o40m",    "o40m", "o8t10m",  "o8t10m",  "u10m",    "u10m",    "u8m", "u8m")
   
  
  
  
   # remove the lines with 0 in weight because due to some landings but during other years we did not keep here.
   land.agg <- land.agg[land.agg$landings !=0,]
  
   landings <- land.agg 
  
  ##---------------------------------------------------------------------
  ## AGGREGATE LANDINGS PER SQUARE---------------------------------------
  ## AND MERGE WITH CODE_SQUARE GIVEN POP SPECIFIC DISTRIBUTION----------
  ##---------------------------------------------------------------------
  
  
   # find out areas
   landings[,c('SI_LATI','SI_LONG')] <- ICESrectangle2LonLat(landings$rectangle, midpoint = TRUE)
   landings$x <- as.numeric(as.character((landings$SI_LONG) ))
   landings$y <-  as.numeric(as.character((landings$SI_LATI) )) 
   landings$code_area   <- longlat_to_ICESareas(landings)
   landings[is.na(landings$code_area) | !(landings$code_area %in% c('22', '23', '24', '25', '26', '27', '28-1', '28-2', '29', '30', '31', '32',
                         'IIIan', 'IIIas', 'IVa', 'IVb', 'IVc') ), 'code_area'] <- 'out' # if out of range,remove...
   landings[(landings$code_area %in% c('IVa', 'IVb', 'IVc')), 'code_area'] <- 'nsea'
   landings[(landings$code_area %in% c('22', '23', '24')), 'code_area'] <- '2224'
   landings[(landings$code_area %in% c('25', '26', '27', '28-1', '28-2', '29','30', '31', '32')), 'code_area'] <- '2532'
   
   # temporary species names
   landings$species   <- paste(landings$species,".", landings$code_area, sep='')
 
  # correct names for special cases (those across management areas)
   landings[ landings$species %in% c("COD.IIIan"), "species"] <- 'COD.nsea'
   landings[ landings$species %in% c("COD.IIIas"), "species"] <- 'COD.kat'
   landings[ landings$species %in% c("HAD.IIIan"), "species"] <- 'HAD.nsea'
   landings[ landings$species %in% c("HER.IIIan", "HER.IIIas", "HER.2224"), "species"] <- 'HER.3a22'
   landings[ landings$species %in% c("SPR.2224", "SPR.2532"), "species"] <- 'SPR.2232'
   landings[ landings$species %in% c("PLE.2224", "PLE.2532"), "species"] <- 'PLE.2232'
   landings[ landings$species %in% c("FLE.2224", "FLE.2532"), "species"] <- 'FLE.2232'
   landings[ landings$species %in% c("TUR.2224", "TUR.2532"), "species"] <- 'TUR.2232'
   landings[ landings$species %in% c("DAB.2224", "DAB.2532"), "species"] <- 'DAB.2232'
   landings[grep('IIIa',  landings$species), "species"] <- # for all other species, correct kask 
     paste( substr(landings[grep('IIIa',  landings$species),'species'], 1,3), ".kask", sep='')

 
 
  # CAUTION !! CAUTION !! CAUTION !! Remove the landings from SIMULATED vessels 
  # unfortunately cannot be  done exactly rigorously because 
  # mismatch between VMS-equipped are not (< and >15m) and DCF segments (<10, 10-12,12-18, etc.)
  
  # first, take a look...
  #unique(landings[landings$country %in% c('DEU', 'DNK'),"vsize"])
  if(general$case_study=="balticonly") unique(landings[landings$country %in% c( 'DNK', 'DEU', 'SWE') ,"vsize"])
  if(general$case_study=="myfish")     unique(landings[landings$country %in% c( 'DNK') ,"vsize"])

  # then remove simulated vessels:
  #idx_to_remove <- which(landings$country %in% c('DEU', 'DNK')  & 
   if(general$case_study=="balticonly"){
     idx_to_remove <- which(landings$country %in% c( 'DNK', 'DEU', 'SWE')  & 
                      landings$vsize %in% c('o18t24m','o24t40m','o40m', 'o15m', 'o12t18m'))
     }
   if(general$case_study=="myfish"){
     idx_to_remove <- which(landings$country %in% c( 'DNK')  & 
                      landings$vsize %in% c('o18t24m','o24t40m','o40m', 'o15m', 'o12t18m'))
     }
  idx_to_keep   <- (1:nrow(landings))[-idx_to_remove]
  
 
 
 
  # here compute the respective part of landings which is
  # from simulated vessels vs. other landings (those recorded into stecf that may be different than the ICES ones...pfff)
  landings_not_simulated_vessels  <- landings[idx_to_keep, ]
  dd1                              <- tapply(landings_not_simulated_vessels$landings, list(landings_not_simulated_vessels$species), sum, na.rm=TRUE)
  landings_simulated_vessels       <- landings[idx_to_remove, ]
  dd2                              <- tapply(landings_simulated_vessels$landings, list(landings_simulated_vessels$species), sum, na.rm=TRUE)
  sp                               <- names(dd2) [names(dd2) %in% names(dd1)]
  compare                          <- cbind.data.frame(not_simu=dd1[names(dd1) %in% sp], simu=dd2[names(dd2) %in% sp]  )
  an                               <- function(x) as.numeric(as.character(x))
  compare$percent_in_simu          <- round(an(compare$simu) /(an(compare$simu)+an(compare$not_simu) )*100 )
  
  ##
  compare$tot <- an(compare$simu)+an(compare$not_simu) # gives the TAC to be put in the pa file for stock parameters step_G_matrix.R  
  ##
  
  pop_names                        <- read.table(file.path(general$main.path,  "popsspe",paste("pop_names_",general$case_study,".txt", sep="")))
  compare[as.character(pop_names[,1]),]
  dd        <- cbind.data.frame(as.character(pop_names[,2]), compare[as.character(pop_names[,1]), 'percent_in_simu' ])
  names(dd) <- c("stock","percent_in_simu")
  dd        <- replace(dd, is.na(dd), 100)
  write.table(dd,    # ordered by species pop_names... 
               file= file.path(general$main.path, "popsspe", 
                  paste('percent_landings_from_simulated_vessels.dat', sep='')), 
                 row.names=FALSE, col.names=TRUE, quote=FALSE)
   
  
  
  
  
  
  # ...and remove
  landings      <- landings[idx_to_keep, ]





  # aggregate
  nm                 <- names(landings)
  idx.col            <- grep('landings', nm) # index columns with species weight
  library(data.table)
  DT                 <- data.table(landings) # library data.table for fast grouping replacing aggregate()
  # AGGREGATE WEIGHT PER SPECIES
  eq1                <- c.listquote( paste ("sum(",nm[idx.col],",na.rm=TRUE)",sep="") )
  DT$landings        <- as.double(DT$landings)
  DT$rectangle       <- as.factor(DT$rectangle)
  DT$species         <- as.factor(DT$species)
  land.agg           <- DT[,eval(eq1),by=list(rectangle, species)]
  land.agg           <- data.frame( land.agg)
  colnames(land.agg) <- c("rectangle", "species", nm[idx.col])


  
  # keep pops of interest only
  pop_names                    <- read.table(file.path(general$main.path,  "popsspe",paste("pop_names_",general$case_study,".txt", sep="")))
  land.agg                     <- land.agg[land.agg$species %in% pop_names[,1],] # keep simulated stocks only 
  land.agg$species             <- factor(land.agg$species)
  levels(land.agg$species)     <- pop_names[,2][ match(levels(land.agg$species), as.character(pop_names[,1]))] # map the name to integer
  land.agg$species             <- as.numeric(as.character(land.agg$species))
  library(doBy)
  land.agg                     <- orderBy(~species, data=land.agg) # library(doBy) # order from 0 to nbstock
 
  
   # STECF data are annual while we need semester data at least 
   # (because different spatail avai key in function of semester)
   #=> assume a division by 2
   land.agg$landings <- land.agg$landings/2  
   land.agg$semester <- 1
   land.agg2 <- land.agg
   land.agg2$semester <- 2
   land.agg <- rbind.data.frame(land.agg, land.agg2)
   
   
   # STECF data are in TONS 
   # while parameterisation for c++ need kilos...
   land.agg$landings <- land.agg$landings *1000 
  
  
  
   # caution myfish stecf fleet sce
   fleetsce <-  data.frame(sce=1:4, namesce=c('stecf_baseline', 'tor1aOr2a', 'tor1bOr2b', 'tor1cOr2c'))
   
    write.table(fleetsce, quote=FALSE,
                 file=file.path(general$main.path.ibm, paste("multiplier_for_fleetsce", general$application,".dat",sep='')), append=FALSE,
                   row.names=FALSE, col.names=TRUE)

  
   # add coordinates
   land.agg <- cbind.data.frame(land.agg, ICESrectangle2LonLat(land.agg$rectangle, midpoint = TRUE))
  
  
    # merge (pop by pop)
    # knowing pop-spe and semester-spe list of nodes
    for(sce in fleetsce$sce){
    for(pop in  1:nrow(pop_names)){
     for(month in c(1:12)){
      if(month %in% c(1:6)) semester <- 1
      if(month %in% c(7:12)) semester <- 2
      
      cat(paste("------------", "\n"))
      cat(paste(pop_names[pop,1], "\n"))
      cat(paste("semester ", semester, "\n"))
        # find out the list of nodes from avai (giving the distribtion of the pop over nodes)
        from_avai_this_semester <- 
                   read.table(file.path("C:","Users","fbas","Documents","GitHub","DISPLACE_input_myfish","popsspe_myfish","static_avai",
                      paste(pop_names[pop,2], 'spe_avai_szgroup_nodes_semester', semester, ".dat", sep='')), header=TRUE)
        lst_nodes_from_avai_this_semester <- unique(from_avai_this_semester[,1])
        code_square_pop_semester <- 
              code_square_per_node[code_square_per_node$pt_graph %in% unique(lst_nodes_from_avai_this_semester),]
        
        # subset this pop this semester
        landings_this_pop <- land.agg[ land.agg$species==pop_names[pop,2] & land.agg$semester==semester, ]
      
        landings_in_pol_sce <- 1 # init (!=0 to avoid 0/)
        
        # find out the  the list of nodes  accepted by the sce
        shapePath <- file.path("C:","Users","fbas","Documents","GitHub","DISPLACE_STECF","Shapes") 
        if(fleetsce$namesce[sce]=="tor1aOr2a" && month %in% c(2, 3)){ # Feb & Mar
           pol                 <- readShapePoly(file.path(shapePath, "1_222324"))
           # test for graph to remove closed areas
           coords              <- SpatialPoints(data.frame(SI_LONG=code_square_pop_semester$SI_LONG,SI_LATI=code_square_pop_semester$SI_LATI))
           idx                 <- over(coords,as(pol,"SpatialPolygons"))
           code_square_pop_semester$idx     <- idx
           code_square_pop_semester         <- code_square_pop_semester[is.na(code_square_pop_semester$idx),]  # keep not included only to dispatch stecf landings on it
           
           # test for stecf landings to get landings in closed areas
           if(nrow(landings_this_pop)!=0) {
              coords              <- SpatialPoints(data.frame(SI_LONG=landings_this_pop$SI_LONG,SI_LATI=landings_this_pop$SI_LATI))
              idx                 <- over(coords,as(pol,"SpatialPolygons"))
              landings_this_pop$idx        <- idx
              landings_in_pol_sce    <- sum(landings_this_pop[!is.na(landings_this_pop$idx), "landings"]) 
             }   else{ landings_in_pol_sce<-1}
           }
        if(fleetsce$namesce[sce]=="tor1bOr2b" && month %in% c(7, 8)){  #Jul & Aug
           pol                 <- readShapePoly(file.path(shapePath, "2_2526"))
           pol2                 <- readShapePoly(file.path(shapePath, "4_24"))
           # test for graph to remove closed areas
           coords              <- SpatialPoints(data.frame(SI_LONG=code_square_pop_semester$SI_LONG,SI_LATI=code_square_pop_semester$SI_LATI))
           idx                 <- over(coords,as(pol,"SpatialPolygons"))
           code_square_pop_semester$idx     <- idx
           idx2                 <- over(coords,as(pol2,"SpatialPolygons"))
           code_square_pop_semester$idx2     <- idx2
           code_square_pop_semester         <- code_square_pop_semester[is.na(code_square_pop_semester$idx) & is.na(code_square_pop_semester$idx2),]  # keep not included only to dispatch stecf landings on it
           
           # test for stecf landings to get landings in closed areas
           if(nrow(landings_this_pop)!=0) {
              coords              <- SpatialPoints(data.frame(SI_LONG=landings_this_pop$SI_LONG,SI_LATI=landings_this_pop$SI_LATI))
              idx                 <- over(coords,as(pol,"SpatialPolygons"))
              landings_this_pop$idx        <- idx
              idx2                 <- over(coords,as(pol2,"SpatialPolygons"))
              landings_this_pop$idx2        <- idx2
              landings_in_pol_sce    <- sum(landings_this_pop[!is.na(landings_this_pop$idx) & !is.na(landings_this_pop$idx2), "landings"]) 
              }  else{ landings_in_pol_sce<-1}
           }
        if(fleetsce$namesce[sce]=="tor1cOr2c" && month %in% c(7, 8)){  #Jul & Aug
           pol                 <- readShapePoly(file.path(shapePath, "3_2728"))
           # test for graph to remove closed areas
           coords              <- SpatialPoints(data.frame(SI_LONG=code_square_pop_semester$SI_LONG,SI_LATI=code_square_pop_semester$SI_LATI))
           idx                 <- over(coords,as(pol,"SpatialPolygons"))
           code_square_pop_semester$idx     <- idx
           code_square_pop_semester         <- code_square_pop_semester[is.na(code_square_pop_semester$idx),]  # keep not included only to dispatch stecf landings on it
           
           # test for stecf landings to get landings in closed areas
           if(nrow(landings_this_pop)!=0) {
              coords              <- SpatialPoints(data.frame(SI_LONG=landings_this_pop$SI_LONG,SI_LATI=landings_this_pop$SI_LATI))
              idx                 <- over(coords,as(pol,"SpatialPolygons"))
              landings_this_pop$idx        <- idx
              landings_in_pol_sce    <- sum(landings_this_pop[!is.na(landings_this_pop$idx), "landings"]) 
              }  else{ landings_in_pol_sce<-1}
           }
           
           
           
        # debug
        if(is.na(landings_in_pol_sce) || landings_in_pol_sce==0)  landings_in_pol_sce <- 1
        
        
        
        # a check 
        map(xlim=c(-5,15),ylim=c(50,60))
        points(code_square_pop_semester[,1:2],col=2)
        
        # find out the number of graph nodes per square (presence for this species this semester) 
        # to later make the dispaching of landings between points
        nb_pt_graph_per_square                          <- tapply(code_square_pop_semester$pt_graph, code_square_pop_semester$code_square, length)
        code_square_pop_semester$code_square            <- factor(code_square_pop_semester$code_square)
        code_square_pop_semester$nb_pt_graph_per_square <- nb_pt_graph_per_square[as.character(code_square_pop_semester$code_square)]
    
        
        
        if(nrow(landings_this_pop)!=0){
         
              
           # do the merging
           merged_this_species_this_quarter <- 
              merge(code_square_pop_semester, landings_this_pop, by.x='code_square', by.y='rectangle')

           # divide the landings evenly between nodes
           merged_this_species_this_quarter$landings <- 
               (merged_this_species_this_quarter$landings / merged_this_species_this_quarter$nb_pt_graph_per_square)  
        
           # add the landings from closed polygons
           merged_this_species_this_quarter$landings <-  merged_this_species_this_quarter$landings + (landings_in_pol_sce / length(merged_this_species_this_quarter$landings)) 
         
           # check if merging did not lost some landings (e.g. if potential landings in square that are not in avai?!..)
           cat(paste("success rate for this pop this semester",
             round(sum(merged_this_species_this_quarter$landings)    /  sum(landings_this_pop$landings) *100,2),
              "% \n"))
        
           # convert landings by semester into landings by month assuming equal distribution 
           merged_this_species_this_quarter$landings <- merged_this_species_this_quarter$landings /6
        
           if(nrow(merged_this_species_this_quarter)==0){
              cat(paste(" no match in the data for this sp this semester  ","\n"))
     
               # ...then, put 0s
               merged_this_species_this_quarter <- cbind(code_square_pop_semester$pt_graph, 0)
               colnames(merged_this_species_this_quarter)  <- c('pt_graph', 'landings')
              }
        
        } else {
        cat(paste(" no other landings in the data for this sp this semester  ","\n"))
     
        # ...then, put 0s
        merged_this_species_this_quarter <- cbind(code_square_pop_semester$pt_graph, 0)
        colnames(merged_this_species_this_quarter)  <- c('pt_graph', 'landings')
        }
  
 
            # save for a multimap in c++  pt_graph / landings kg per month (to be used every months)
            # to fill in the pop attribute
        write.table(round(merged_this_species_this_quarter[, c('pt_graph', 'landings')]),
               file= file.path(general$main.path, "popsspe", 
                  paste(pop_names[pop,2], 'spe_stecf_oth_land_per_month_per_node_month', month, "_fleetsce",sce,".dat", sep='')), 
                 row.names=FALSE, col.names=TRUE, quote=FALSE)
   
         
  }}
  } # end fleet sce
  
  # check 
  cc <-read.table(file.path(general$main.path, "popsspe", "10spe_stecf_oth_land_per_month_per_node_month1_fleetsce2.dat"), header=TRUE)
  dd <-read.table(file.path(general$main.path, "popsspe", "10spe_stecf_oth_land_per_month_per_node_month2_fleetsce2.dat"), header=TRUE)