store corrected lat long for knots

.Rbuildignore

@@ -1,4 +1,6 @@
 ^\.git$
 ^\.travis\.yml$
 ^shiny$
-^examples$
+^examples$
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -3,4 +3,6 @@
 *.so
 shiny/database/
 shiny/testing/
-shiny/global_coverage.png
+shiny/global_coverage.png
+.Rproj.user
+*.Rproj

.travis.yml

@@ -45,6 +45,43 @@ r_packages:
   - mapdata
   - TMB
 
+cache: packages
+
+sudo: required
+
+dist: trusty
+
+addons:
+  postgresql: "9.6"
+
+before_install:
+  - sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable --yes
+  - sudo apt-get --yes --force-yes update -qq
+  - sudo apt-get install --yes libudunits2-dev libproj-dev libgeos-dev libgdal-dev
+  - # install postgis from source:
+  - sudo apt-get --yes install libjson-c-dev postgresql-server-dev-9.6
+  - wget http://download.osgeo.org/postgis/source/postgis-2.3.2.tar.gz
+  - (mv postgis* /tmp; cd /tmp; tar xzf postgis-2.3.2.tar.gz)
+  - (cd /tmp/postgis-2.3.2 ; ./configure; make; sudo make install)
+  - # activate liblwgeom:
+  - sudo ldconfig
+  - # create postgis databases:
+  - sudo service postgresql restart
+  - createdb postgis
+  - psql -d postgis -c "CREATE EXTENSION postgis;"
+  - psql -d postgis -c "GRANT CREATE ON DATABASE postgis TO travis"
+  - createdb empty
+  - psql -d empty -c "CREATE EXTENSION postgis;"
+
+warnings_are_errors: true
+
+after_success:
+  - dropdb postgis
+  - createdb postgis
+  - psql -d postgis -c "CREATE EXTENSION postgis;"
+  - psql -d postgis -c "GRANT CREATE ON DATABASE postgis TO travis"
+  - Rscript -e 'covr::codecov()'
+
 r_github_packages:
   - kaskr/TMB_contrib_R/TMBhelper
 

DESCRIPTION

@@ -12,7 +12,7 @@ Description: SpatialDeltaGLMM is an R package for implementing a spatial delta-
     functions and model-comparison tools, and (3) is intended to improve analysis
     speed, replicability, peer-review, and interpretation of index standardization.
     methods
-Imports: 
+Imports:
     graphics,
     utils,
     mapproj,
@@ -31,17 +31,21 @@ Imports:
     maps,
     mapdata,
     TMB,
-    MatrixModels  
-Depends: 
+    MatrixModels  ,
+    dplyr,
+    purrr,
+    sf,
+    rgdal
+Depends:
     R (>= 3.1.0),
-Suggests: 
+Suggests:
     INLA,
     testthat
-Additional_repositories: 
+Additional_repositories:
     https://www.math.ntnu.no/inla/R/stable
 License: GPL-2
 LazyData: yes
 BuildVignettes: yes
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1
 URL: http://github.com/nwfsc-assess/geostatistical_delta-GLMM
 BugReports: http://github.com/nwfsc-assess/geostatistical_delta-GLMM/issues

NAMESPACE

@@ -36,6 +36,7 @@ export(Spatial_Information_Fn)
 export(Summarize)
 export(Vessel_Fn)
 export(bilinear_interp)
+export(convert_utm_to_ll_fn)
 export(match_strata_fn)
 export(plot_residuals)
 export(summary_nwfsc)

R/Spatial_Information_Fn.R

@@ -1,4 +1,6 @@
 
+
+
 #' Build objects related to spatial information
 #'
 #' \code{Spatial_Information_Fn} builds a tagged list with all the spatial information needed for \code{Data_Fn}
@@ -17,78 +19,174 @@
 #'   \item{MeshList}{A tagged list with inputs related to the SPDE mesh}
 #'   \item{GridList}{A tagged list with inputs related to the 2D AR1 grid}
 #'   \item{a_xl}{A data frame with areas for each knot and each strattum}
-#'   \item{loc_UTM}{A data frame with the converted UTM coordinates for each sample}
+#'   \item{loc_x}{The UTM location for each knot}
+#'   \item{loc_x_lat_long}{A data frame with the corrected lat long for each knot}
 #'   \item{Kmeans}{Output from \code{Calc_Kmeans} with knots for a triangulated mesh}
 #'   \item{knot_i}{The knot associated with each sample}
 #'   \item{Method}{The Method input (for archival purposes)}
-#'   \item{loc_x}{The UTM location for each knot}
 #' }
 
 #' @export
-Spatial_Information_Fn = function( n_x, Lon, Lat, Extrapolation_List, Method="Mesh", grid_size_km=50, grid_size_LL=1, ... ){
-
+Spatial_Information_Fn = function(n_x,
+                                  Lon,
+                                  Lat,
+                                  Extrapolation_List,
+                                  Method = "Mesh",
+                                  grid_size_km = 50,
+                                  grid_size_LL = 1,
+                                  ...) {
   # Convert to an Eastings-Northings coordinate system
-  if( Method=="Spherical_mesh" ){
-    loc_i = data.frame( 'Lon'=Lon, 'Lat'=Lat )
+  if (Method == "Spherical_mesh") {
+    loc_i = data.frame('Lon' = Lon, 'Lat' = Lat)
     # Bounds for 2D AR1 grid
-    Grid_bounds = (grid_size_km/110) * apply(Extrapolation_List$Data_Extrap[,c('Lon','Lat')]/(grid_size_km/110), MARGIN=2, FUN=function(vec){trunc(range(vec))+c(0,1)})
+    Grid_bounds = (grid_size_km / 110) * apply(
+      Extrapolation_List$Data_Extrap[, c('Lon', 'Lat')] / (grid_size_km / 110),
+      MARGIN = 2,
+      FUN = function(vec) {
+        trunc(range(vec)) + c(0, 1)
+      }
+    )
 
     # Calculate k-means centroids
-    Kmeans = SpatialDeltaGLMM::Calc_Kmeans(n_x=n_x, loc_orig=loc_i[,c("Lon", "Lat")], ... )
+    Kmeans = SpatialDeltaGLMM::Calc_Kmeans(n_x = n_x, loc_orig = loc_i[, c("Lon", "Lat")], ...)
 
     # Calculate grid for 2D AR1 process
-    loc_grid = expand.grid( 'Lon'=seq(Grid_bounds[1,1],Grid_bounds[2,1],by=grid_size_LL), 'Lat'=seq(Grid_bounds[1,2],Grid_bounds[2,2],by=grid_size_LL) )
-    Which = sort(unique(RANN::nn2(data=loc_grid, query=Extrapolation_List$Data_Extrap[which(Extrapolation_List$Area_km2_x>0),c('Lon','Lat')], k=1)$nn.idx[,1]))
+    loc_grid = expand.grid(
+      'Lon' = seq(Grid_bounds[1, 1], Grid_bounds[2, 1], by = grid_size_LL),
+      'Lat' = seq(Grid_bounds[1, 2], Grid_bounds[2, 2], by = grid_size_LL)
+    )
+    Which = sort(unique(
+      RANN::nn2(
+        data = loc_grid,
+        query = Extrapolation_List$Data_Extrap[which(Extrapolation_List$Area_km2_x >
+                                                       0), c('Lon', 'Lat')],
+        k = 1
+      )$nn.idx[, 1]
+    ))
     loc_grid = loc_grid[Which,]
-    grid_num = RANN::nn2( data=loc_grid, query=loc_i, k=1)$nn.idx[,1]
+    grid_num = RANN::nn2(data = loc_grid,
+                         query = loc_i,
+                         k = 1)$nn.idx[, 1]
   }
-  if( Method %in% c("Mesh","Grid") ){
-    loc_i = SpatialDeltaGLMM::Convert_LL_to_UTM_Fn( Lon=Lon, Lat=Lat, zone=Extrapolation_List$zone, flip_around_dateline=Extrapolation_List$flip_around_dateline )                                                         #$
-    loc_i = cbind( 'E_km'=loc_i[,'X'], 'N_km'=loc_i[,'Y'])
+
+
+  if (Method %in% c("Mesh", "Grid")) {
+    loc_i = SpatialDeltaGLMM::Convert_LL_to_UTM_Fn(
+      Lon = Lon,
+      Lat = Lat,
+      zone = Extrapolation_List$zone,
+      flip_around_dateline = Extrapolation_List$flip_around_dateline
+    )                                                         #$
+    loc_i = cbind('E_km' = loc_i[, 'X'], 'N_km' = loc_i[, 'Y'])
     # Bounds for 2D AR1 grid
-    Grid_bounds = grid_size_km * apply(Extrapolation_List$Data_Extrap[,c('E_km','N_km')]/grid_size_km, MARGIN=2, FUN=function(vec){trunc(range(vec))+c(0,1)})
+    Grid_bounds = grid_size_km * apply(
+      Extrapolation_List$Data_Extrap[, c('E_km', 'N_km')] / grid_size_km,
+      MARGIN = 2,
+      FUN = function(vec) {
+        trunc(range(vec)) + c(0, 1)
+      }
+    )
 
     # Calculate k-means centroids
-    Kmeans = SpatialDeltaGLMM::Calc_Kmeans(n_x=n_x, loc_orig=loc_i[,c("E_km", "N_km")], ... )
+    Kmeans = SpatialDeltaGLMM::Calc_Kmeans(n_x = n_x, loc_orig = loc_i[, c("E_km", "N_km")], ...)
 
     # Calculate grid for 2D AR1 process
-    loc_grid = expand.grid( 'E_km'=seq(Grid_bounds[1,1],Grid_bounds[2,1],by=grid_size_km), 'N_km'=seq(Grid_bounds[1,2],Grid_bounds[2,2],by=grid_size_km) )
-    Which = sort(unique(RANN::nn2(data=loc_grid, query=Extrapolation_List$Data_Extrap[which(Extrapolation_List$Area_km2_x>0),c('E_km','N_km')], k=1)$nn.idx[,1]))
+    loc_grid = expand.grid(
+      'E_km' = seq(Grid_bounds[1, 1], Grid_bounds[2, 1], by = grid_size_km),
+      'N_km' = seq(Grid_bounds[1, 2], Grid_bounds[2, 2], by = grid_size_km)
+    )
+    Which = sort(unique(
+      RANN::nn2(
+        data = loc_grid,
+        query = Extrapolation_List$Data_Extrap[which(Extrapolation_List$Area_km2_x >
+                                                       0), c('E_km', 'N_km')],
+        k = 1
+      )$nn.idx[, 1]
+    ))
     loc_grid = loc_grid[Which,]
-    grid_num = RANN::nn2( data=loc_grid, query=loc_i, k=1)$nn.idx[,1]
+    grid_num = RANN::nn2(data = loc_grid,
+                         query = loc_i,
+                         k = 1)$nn.idx[, 1]
   }
 
   # Calc design matrix and areas
-  if( Method=="Grid" ){
+  if (Method == "Grid") {
     knot_i = grid_num
     loc_x = loc_grid
   }
-  if( Method %in% c("Mesh","Spherical_mesh") ){
+  if (Method %in% c("Mesh", "Spherical_mesh")) {
     knot_i = Kmeans$cluster
     loc_x = Kmeans$centers
   }
-  PolygonList = Calc_Polygon_Areas_and_Polygons_Fn( loc_x=loc_x, Data_Extrap=Extrapolation_List[["Data_Extrap"]], a_el=Extrapolation_List[["a_el"]])
-  a_xl = PolygonList[["a_xl"]]
 
+  eastings <- loc_x[, 'E_km']
+
+  northings <- loc_x[, 'N_km']
+
+  zone <- Extrapolation_List$zone
+
+  flip_around_dateline <- Extrapolation_List$flip_around_dateline
+  loc_x_lat_long <-
+    convert_utm_to_ll_fn(
+      eastings = eastings,
+      northings = northings,
+      zone = zone,
+      flip_around_dateline = flip_around_dateline
+    )
+
+  loc_x_lat_long$knot <- 1:nrow(loc_x_lat_long)
+
+  if (flip_around_dateline == T){
+    warning("knots have been flipped around dateline; eastings locations and UTM zone are transformed - see loc_x_lat_long for corrected locations")
+  }
+
+  PolygonList = Calc_Polygon_Areas_and_Polygons_Fn(loc_x = loc_x,
+                                                   Data_Extrap = Extrapolation_List[["Data_Extrap"]],
+                                                   a_el = Extrapolation_List[["a_el"]])
+  a_xl = PolygonList[["a_xl"]]
   # Convert loc_x back to location in lat-long coordinates loc_x_LL
   # if zone=NA or NULL, then it automatically detects appropriate zone
-  #tmpUTM = cbind('PID'=1,'POS'=1:nrow(loc_x),'X'=loc_x[,'E_km'],'Y'=loc_x[,'N_km'])
-  #attr(tmpUTM,"projection") = "UTM"
-  #attr(tmpUTM,"zone") = Extrapolation_List$zone
-  #loc_x_LL = PBSmapping::convUL(tmpUTM)                                                         #$
+  # tmpUTM = cbind('PID'=1,'POS'=1:nrow(loc_x),'X'=loc_x[,'E_km'],'Y'=loc_x[,'N_km'])
+  # attr(tmpUTM,"projection") = "UTM"
+  # attr(tmpUTM,"zone") = Extrapolation_List$zone
+  # loc_x_LL = PBSmapping::convUL(tmpUTM)                                                         #$
 
   # Make mesh and info for anisotropy  SpatialDeltaGLMM::
-  MeshList = Calc_Anisotropic_Mesh( Method=Method, loc_x=Kmeans$centers )
+  MeshList = Calc_Anisotropic_Mesh(Method = Method, loc_x = Kmeans$centers)
 
   # Make matrices for 2D AR1 process
-  Dist_grid = dist(loc_grid, diag=TRUE, upper=TRUE)
-  M0 = as( ifelse(as.matrix(Dist_grid)==0, 1, 0), "dgTMatrix" )
-  M1 = as( ifelse(as.matrix(Dist_grid)==grid_size_km, 1, 0), "dgTMatrix" )
-  M2 = as( ifelse(as.matrix(Dist_grid)==sqrt(2)*grid_size_km, 1, 0), "dgTMatrix" )
-  if( Method=="Spherical_mesh" ) GridList = list("M0"=M0, "M1"=M1, "M2"=M2, "grid_size_km"=grid_size_LL)
-  if( Method %in% c("Mesh","Grid") ) GridList = list("M0"=M0, "M1"=M1, "M2"=M2, "grid_size_km"=grid_size_km)
+  Dist_grid = dist(loc_grid, diag = TRUE, upper = TRUE)
+  M0 = as(ifelse(as.matrix(Dist_grid) == 0, 1, 0), "dgTMatrix")
+  M1 = as(ifelse(as.matrix(Dist_grid) == grid_size_km, 1, 0), "dgTMatrix")
+  M2 = as(ifelse(as.matrix(Dist_grid) == sqrt(2) * grid_size_km, 1, 0), "dgTMatrix")
+  if (Method == "Spherical_mesh")
+    GridList = list(
+      "M0" = M0,
+      "M1" = M1,
+      "M2" = M2,
+      "grid_size_km" = grid_size_LL
+    )
+  if (Method %in% c("Mesh", "Grid"))
+    GridList = list(
+      "M0" = M0,
+      "M1" = M1,
+      "M2" = M2,
+      "grid_size_km" = grid_size_km
+    )
 
   # Return
-  Return = list("MeshList"=MeshList, "GridList"=GridList, "a_xl"=a_xl, "loc_i"=loc_i, "Kmeans"=Kmeans, "knot_i"=knot_i, "Method"=Method, "loc_x"=loc_x, "PolygonList"=PolygonList, "NN_Extrap"=PolygonList$NN_Extrap)
-  return( Return )
+  Return = list(
+    "MeshList" = MeshList,
+    "GridList" = GridList,
+    "a_xl" = a_xl,
+    "loc_i" = loc_i,
+    "Kmeans" = Kmeans,
+    "knot_i" = knot_i,
+    "Method" = Method,
+    "loc_x" = loc_x,
+    "loc_x_lat_long" = loc_x_lat_long,
+    "PolygonList" = PolygonList,
+    "NN_Extrap" = PolygonList$NN_Extrap
+  )
+  return(Return)
 }