## Warning: package 'tidyr' was built under R version 3.6.2
## Warning: package 'DT' was built under R version 3.6.2
## Warning: package 'anytime' was built under R version 3.6.2
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date

Topographic survey data collected in Ackerson Meadow using two Emlid Reach RTK GPS rovers. Surveying was conducted over two days by J. Shaw and E. Gage

Ackerson-shaw080319.shp

contains all float or fix points (except for one labeled as “none” for solution, which is eliminated here)

raw.d1a.shaw <- st_read("./survey_data_2019/raw/Ackerson-shaw080319.shp")
## Reading layer `Ackerson-shaw080319' from data source `D:\Dropbox\AEI\Projects\AckersonMeadow\survey_data_2019\raw\Ackerson-shaw080319.shp' using driver `ESRI Shapefile'
## Simple feature collection with 270 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -119.8375 ymin: 0 xmax: 0 ymax: 37.84184
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
raw.d1a.shaw <- raw.d1a.shaw %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names()%>% 
  mutate(job = "Shaw") %>% 
  mutate(file = "Ackerson-shaw080319.shp")

raw.d1a.shaw <- 
  raw.d1a.shaw %>% 
  filter(solution_s != "NONE")

raw.d1a.shaw %>%
  as_tibble() %>% 
  janitor::tabyl(solution_s) %>% 
  gt() %>% 
  tab_header(title = "Ackerson-shaw080319.shp")
Ackerson-shaw080319.shp
solution_s n percent
FIX 150 0.5576208
FLOAT 119 0.4423792
NONE 0 0.0000000

Ackerson gage 01.shp

contains all ‘single’ solution points. To be post-processed using base logs

raw.d1a.gage <- st_read("./survey_data_2019/raw/Ackerson gage 01.shp")
## Reading layer `Ackerson gage 01' from data source `D:\Dropbox\AEI\Projects\AckersonMeadow\survey_data_2019\raw\Ackerson gage 01.shp' using driver `ESRI Shapefile'
## Simple feature collection with 117 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -119.8308 ymin: 37.84117 xmax: -119.8269 ymax: 37.84174
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
raw.d1a.gage <- raw.d1a.gage %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names()%>% 
  mutate(job = "Gage") %>% 
  mutate(file = "Ackerson gage 01.shp")


raw.d1a.gage %>%
  as_tibble() %>%
  mutate(dt = lubridate::date(collection)) %>% 
  janitor::tabyl(solution_s, dt) %>% 
  gt() %>% 
  tab_header(title = "Ackerson gage 01.shp")
Ackerson gage 01.shp
solution_s 2019-08-02 2019-08-03
SINGLE 72 45 
mapview(raw.d1a.gage, zcol = 'solution_s')

Ackerson-Shaw.shp contains a mix of float, fix, and some ‘single’ solution points

raw.d1b.shaw <- st_read("./survey_data_2019/raw/Ackerson-Shaw.shp")
## Reading layer `Ackerson-Shaw' from data source `D:\Dropbox\AEI\Projects\AckersonMeadow\survey_data_2019\raw\Ackerson-Shaw.shp' using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -119.8307 ymin: 37.84124 xmax: -119.8269 ymax: 37.84212
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
raw.d1b.shaw <- raw.d1b.shaw %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names()%>% 
  mutate(job = "Shaw") %>% 
  mutate(file = "Ackerson-Shaw.shp")

raw.d1b.shaw %>%
  as_tibble() %>% 
  mutate(dt = lubridate::as_datetime(collection)) %>% 
  mutate(dt_pac = with_tz(dt, tzone = "America/Los_Angeles")) %>%
  janitor::tabyl(solution_s, dt_pac) %>% 
  gt() %>% 
  tab_header(title = "Ackerson-Shaw.shp")
Ackerson-Shaw.shp
solution_s 2019-08-02 16:06:47 2019-08-02 16:07:39 2019-08-02 16:08:43 2019-08-02 16:09:28 2019-08-02 16:10:13 2019-08-02 16:11:50 2019-08-02 16:12:35 2019-08-02 16:13:23 2019-08-02 16:14:22 2019-08-02 16:15:04 2019-08-02 16:15:55 2019-08-02 16:17:12 2019-08-02 16:17:36 2019-08-02 16:18:28 2019-08-02 16:18:51 2019-08-02 16:19:48 2019-08-02 16:20:15 2019-08-02 16:21:13 2019-08-02 16:22:25 2019-08-02 16:23:56 2019-08-02 16:25:22 2019-08-02 16:27:04 2019-08-02 16:28:26 2019-08-02 16:29:08 2019-08-02 16:29:42 2019-08-02 16:33:17 2019-08-02 16:34:21 2019-08-02 16:35:05 2019-08-02 16:40:45 2019-08-02 16:43:35 2019-08-02 16:43:58 2019-08-02 16:44:22 2019-08-02 16:45:01 2019-08-02 16:45:33 2019-08-02 16:46:13 2019-08-02 16:47:06 2019-08-02 16:48:22 2019-08-02 16:48:40 2019-08-02 16:49:41 2019-08-02 16:50:34 2019-08-02 16:51:58 2019-08-02 16:52:20 2019-08-02 16:53:26 2019-08-02 16:53:51 2019-08-02 16:54:21 2019-08-02 16:54:42 2019-08-02 16:55:41 2019-08-02 16:56:53 2019-08-02 16:58:04 2019-08-02 16:58:42 2019-08-02 16:59:15 2019-08-02 16:59:59 2019-08-02 17:00:35 2019-08-02 17:01:26 2019-08-02 17:02:00 2019-08-02 17:02:29 2019-08-02 17:03:01 2019-08-02 17:03:39 2019-08-02 17:04:19 2019-08-02 17:06:43 2019-08-02 17:07:10 2019-08-02 17:07:32 2019-08-02 17:08:11 2019-08-02 17:08:28 2019-08-02 17:09:50 2019-08-02 17:10:51 2019-08-02 17:11:32 2019-08-02 17:12:31 2019-08-02 17:13:02 2019-08-02 17:13:48 2019-08-02 17:14:22 2019-08-02 17:15:11 2019-08-02 17:15:36 2019-08-02 17:16:34 2019-08-02 17:17:22 2019-08-02 17:17:43 2019-08-02 17:18:13 2019-08-02 17:18:27 2019-08-02 17:19:07 2019-08-02 17:20:15 2019-08-02 17:21:09 2019-08-02 17:22:09 2019-08-02 17:24:30 2019-08-02 17:26:19 2019-08-02 17:26:55 2019-08-02 17:29:20 2019-08-02 17:30:42
FIX 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
FLOAT 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 0 1 1 1 0 1 1 1 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1
SINGLE 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 
# Ackerson-Shaw.shp

mapviewOptions(basemaps = c("Esri.WorldImagery", "OpenStreetMap"))
mapview(raw.d1b.shaw, zcol = 'solution_s')
# read in raw survey data
raw.d2.gage <- st_read("./survey_data_2019/raw/Ackerson gage2.shp", quiet = TRUE)

raw.d2.gage <- raw.d2.gage %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names() %>% 
  mutate(job = "Gage") %>% 
  mutate(file = "Ackerson gage2.shp") %>% 
  mutate(antenna_he = 2.165) # change RH to height in field book

raw.d2.shaw <- st_read("./survey_data_2019/raw/Ackerson-Shaw 080419.shp")
## Reading layer `Ackerson-Shaw 080419' from data source `D:\Dropbox\AEI\Projects\AckersonMeadow\survey_data_2019\raw\Ackerson-Shaw 080419.shp' using driver `ESRI Shapefile'
## Simple feature collection with 327 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -119.8461 ymin: 37.83585 xmax: -119.8375 ymax: 37.83956
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
raw.d2.shaw <- raw.d2.shaw %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names()%>% 
  mutate(job = "Shaw") %>% 
  mutate(file = "Ackerson-Shaw 080419.shp")

mapview(raw.d2.gage, zcol = 'solution_s')
mapview(raw.d2.shaw, zcol = 'solution_s')
raw.d3.gage <- st_read("./survey_data_2019/raw/Ackerson gage 3.shp", quiet = TRUE)

raw.d3.gage <- raw.d3.gage %>% 
  mutate(name = tolower(name)) %>%
  janitor::clean_names() %>% 
  mutate(job = "Gage") %>% 
  mutate(antenna_he = 2.165) %>% 
  mutate(file = "Ackerson gage 3.shp")#

raw.d3.gage %>% 
  glimpse()
## Observations: 368
## Variables: 12
## $ antenna_he <dbl> 2.165, 2.165, 2.165, 2.165, 2.165, 2.165, 2.165, 2.165, ...
## $ collection <fct> 2019-08-04 16:56:01 UTC, 2019-08-04 16:57:14 UTC, 2019-0...
## $ collecti_1 <fct> 2019-08-04 16:55:54 UTC, 2019-08-04 16:57:11 UTC, 2019-0...
## $ lateral_rm <dbl> 0.0167, 0.0049, 0.0053, 0.0034, 0.0378, 0.0164, 0.0420, ...
## $ name       <chr> "bm usda", "bm nps", "bm nps", "bm usda", "tb", "tb", "t...
## $ projected  <fct> -119.831621116 37.8405594082 1393.8686973, -119.83162279...
## $ rms        <fct> 0.0039582103794 0.0162589911839 0.00230440200514, 0.0044...
## $ sample_cou <dbl> 37, 15, 10, 8, 10, 8, 9, 8, 8, 8, 10, 7, 9, 8, 18, 7, 7,...
## $ solution_s <fct> FIX, FIX, FIX, FIX, FIX, FIX, FIX, FIX, FIX, FIX, FIX, F...
## $ geometry   <POINT [°]> POINT (-119.8316 37.84056), POINT (-119.8316 37.84...
## $ job        <chr> "Gage", "Gage", "Gage", "Gage", "Gage", "Gage", "Gage", ...
## $ file       <chr> "Ackerson gage 3.shp", "Ackerson gage 3.shp", "Ackerson ...

merge points

### bind d1
survey.d1.raw <- rbind(raw.d1a.gage, raw.d1a.shaw)
survey.d1.raw <- rbind(survey.d1.raw, raw.d1b.shaw)

### bind d2
survey.d2.raw <- rbind(raw.d2.gage,raw.d2.shaw)


### bind d1 and d2
survey.d1d2.raw <- rbind(survey.d1.raw, survey.d2.raw)

survey.d1d2d3.raw <- rbind(survey.d1d2.raw,raw.d3.gage)


mapview(survey.d1d2.raw)
mapview(survey.d1d2d3.raw, zcol = 'file')
## separate projected coordinates, rmse values
survey.d1d2d3.raw <- survey.d1d2d3.raw %>% 
  separate(projected, c("long", "lat", "elev_m_raw"), sep = " ") %>%
  separate(rms, c("rms_x", "rms_y", "rms_z"), sep = " ")

## convert char to numeric, calculate the adjusted elev
survey.d1d2d3.raw <- survey.d1d2d3.raw %>%
  mutate_at(.,.vars = c("long", "lat", "elev_m_raw", "rms_x", "rms_y", "rms_z"),.funs = as.numeric) %>% 
  mutate(elev_m_adj = elev_m_raw - antenna_he) %>% 
  mutate(dt = as_datetime(collection)) 


mapview(survey.d1d2d3.raw, zcol = 'elev_m_adj')
survey.d1d2d3.raw %>% 
  # glimpse()
  ggplot(aes(dt,antenna_he)) +
  geom_point(aes(color = job))

survey.d1d2d3.raw %>% 
  mutate(dt = as_datetime(collection)) %>% 
  ggplot(aes(dt,antenna_he)) +
  geom_point(aes(color = file))

survey.d1d2d3.raw %>% 
  ggplot(aes(x = elev_m_adj)) +
  geom_density()

survey.d1d2d3.raw %>% 
  ggplot(aes(x = file, y = elev_m_adj)) +
  geom_boxplot() 

plt1 <- survey.d1d2d3.raw %>% 
  ggplot(aes(x = file, y = elev_m_adj)) +
  geom_boxplot() +
  geom_jitter(alpha=.2)
plt1

plotly::ggplotly(plt1)
plt2 <- survey.d1d2d3.raw %>% 
  ggplot(aes(x = dt, y = elev_m_adj)) +
  geom_point(aes(color = file),alpha=.2) +
  facet_wrap(~job, ncol=1)
plt2

plotly::ggplotly(plt2)

Count of points by solution type

## 
survey.d1d2d3.raw %>%
  # glimpse()
  janitor::tabyl(solution_s) %>% 
  gt()
solution_sgeometry n percent
SINGLE 150 0.08886256
FIX 682 0.40402844
FLOAT 856 0.50710900
NONE 0 0.00000000 
# survey.d1d2.fixfloat
survey.fixfloat <- survey.d1d2d3.raw %>% 
  filter(solution_s != "SINGLE")

Map: day 1 and day 2 raw points

mapviewOptions(basemaps = c("Esri.WorldImagery", "OpenStreetMap"))
mapview(survey.d1d2.raw, zcol = "job")

Map: day 1 and day 2 raw solution

mapview(survey.d2.raw, zcol = "solution_s")
mapviewOptions(basemaps = c("Esri.WorldImagery", "OpenStreetMap"))

Map: Fix/float points

survey.fixfloat %>% 
  # filter(rms_z > .07) %>% 
  mapview(., zcol = "elev_m_adj")
mapviewOptions(basemaps = c("Esri.WorldImagery", "OpenStreetMap"))

Map: High RMS Z (>0.07)

survey.fixfloat %>% 
  filter(rms_z > .07) %>% 
  mapview(., zcol = "rms_z")
mapviewOptions(basemaps = c("Esri.WorldImagery", "OpenStreetMap"))
pl1 <- survey.d2.raw %>% 
  ggplot(aes(rms_x, rms_y)) +
  geom_point(aes(color = rms_z), shape = 16 ) +
  scale_color_viridis(option = 'B') +
  geom_text(data=subset(survey.raw, rms_z > 0.11 | rms_x > 0.025),
            aes(rms_x, rms_y, label=name))
library(plotly)
ggplotly(pl1)
survey.d1d2d3.raw %>%
  # names()
  ggplot(aes(rms_x, rms_y)) +
  geom_point(aes(fill = rms_z), shape = 21 ) +
  scale_fill_viridis(option = 'B') #+

  # ggrepel::geom_text_repel(data=subset(survey.d1d2.raw, rms_z > 0.1 & rms_x > 0.025),
  #           aes(rms_x, rms_y, label=name))
st_write(survey.d1d2.fixfloat, "./Ackerson_d1d2_fixfloat.shp")