This and other labs in this course are based on those developed by C.J. Gabbe and Iris Stewart-Frey.

Today we are going to be dealing data cleaning and what is known as data wrangling within data tables, we will also convert the information to a spatial dataset. Later on this quarter we will use a great package called tidycensus to apply similar techniques to spatial data right from the start

Steps in today’s lab

Read csv and explore
Convert csv to sf
Project sf data
Remove duplicates
Filter observations
Remove unnecessary columns
Rename variables
Add new variables (mutate)
Manage missing values appropriately
Re-write steps 4-9 using pipes and efficient coding
Organize your RScript and save

Prep work

Create a ‘Module3_Lab’ folder in the drive of your choice (preferably somewhere that’s backed up into the cloud)

Create a new RScript by going to File > New File > RScript. Save the RScript to your ‘Module3_Lab’ folder with a name of your choosing. Create a header in your script.

Load the packages we’re going to use today and set your working directory:

library(tidyverse)

## Warning: package 'ggplot2' was built under R version 4.3.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

setwd("~/Google Drive/My Drive/ENVS 117/04_Labs/Module3_Lab")

setwd("~/Desktop/ENVS117")

1. Read csv and explore

Today we’re going to explore data about well contamination in the State of California’s GAMA database. At SCU our water and climate lab has selected all of the available nitrate test results for wells in San Joaquin, Merced, Stanislaus, and Tuolumne Counties, unfornutaley we need to clean some of the data up like removing NA values as well as getting rid of some duplicate data. Download it here, then read it into R:

WELLS_CV_NO3 <- read_csv("GAMA_2021_Well_Data_01.csv")

## New names:
## Rows: 56452 Columns: 16
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (8): well id, chemical id, CHEMICAL, UNIT, data source, WELL_TYPE, DATE,... dbl
## (8): ...1, GM_Results, GM_Longitude, GM_Latitude, MCL, WELL_DEPTH_FT, TO...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`

There are a lot of variables in this dataset! Let’s refer to the Data Dictionary and the Data Description to figure out what some of these fields are all about.

2. Convert csv to sf

You’ll notice that the dataset includes longitude and latitude for each well. This is very helpful, and makes it possible to convert into sf (simple features) format.

What does the st_as_sf command do? What are the options that I’ve included in the code below? This function converts a dataframe to a simple feature using spatial data in a column or columns. It also as an argument called crs which stands for coordinate refrence system. Just like we did in ENVS 116 we are projecting our data. Use infomration from typing ? st_as_sf

WELLS_CV_NO3.sf <- st_as_sf(WELLS_CV_NO3, coords = c("GM_Longitude", "GM_Latitude"), 
                             crs = 4269, agr = "constant") # Convert to sf, crs = NAD 1983

What is EPSG:4269? Google it to find out.

3. Project data

We’re going to project our nitrate well contamination dataset into the California Albers projection, which is the most commonly-used statewide projection in California. Look at this overview or to learn about how projections are represented in R. You can find the EPSG code for California Albers by looking at Spatial References, or by searching for “ESPG code for California Albers R” online, i.e. here. Then, we’ll project our data:

WELLS_CV_NO3.sf <- st_transform(WELLS_CV_NO3.sf, crs = 3310)

#3a Lets try plotting the plot and ggplot function data. Notice again that we have to remove the # from in front of the code

plot(WELLS_CV_NO3.sf)

## Warning: plotting the first 9 out of 14 attributes; use max.plot = 14 to plot
## all

Now lets try and plot one single variable which will be the GM_Results variable

library(ggplot2)
ggplot() + 
  geom_sf(data = WELLS_CV_NO3.sf, aes(fill = GM_Results))  + 
  scale_y_continuous(breaks = 0:100)

We have some values that are exceptionally high so lets subset for all vlues less than 50 mg/l.

WELLS_CV_NO3.sfSubset <-subset(WELLS_CV_NO3.sf, GM_Results <50)

Now lets map it out for values 0 to 100. to do so wwe will use the subset sf feature and add the scale_color_viridis to add a color scale. you can try to set option to C or D as well if you would like a different color ramp

ggplot() + 
  geom_sf(data = WELLS_CV_NO3.sfSubset, aes(color = GM_Results))  + 
  scale_color_viridis_c(option = "B")

Now lets return back to our original sf feature

4. Remove duplicates

Sometimes datasets include duplicate values. These can be across all variables (like someone accidentally added the same observation twice) or they could be for a single variable (two observations have the exact same address). The distinct command lets us remove the duplicates for all observations. What does the option .keep_all = TRUE do? Here you can once again run in the console ? plus the function in this case [distinct]

WELLS_CV_NO3.sf <- distinct(WELLS_CV_NO3.sf, .keep_all = TRUE)

5. Filter observations

We’re most interested in the data from counties that are in the lower, agricultural portions of the San Joaquin Valley, so let’s filter all the wells that are in Merced, and Stanislaus Counties. The filter command lets us select only those observations that fit our criteria:

WELLS_CV_NO3.sf <- dplyr::filter(WELLS_CV_NO3.sf, county == "merced" | county == "stanislaus" )

You’ll see that I used dplyr::filter above. Is there any difference between filter and dplyr::filter? Hint: type ?filter into the console.

6. Remove unnecessary columns

This dataset includes several variables that are not going to be relevant to us, for example, there are two columns that indicate the contaminant (nitrate). So let’s keep only those columns that we need, specifically:

GM_Results
well id
chemical id
UNIT
WELL_DEPTH_FT
data source
WELL_TYPE
TOP_OF_SCREEN
BOTTOM_OF_SCREEN
DATE
county
geometry

We can select variables by column number or by name. I often use column number. You can get a list of variables and their associated column numbers using names, as you can see below. After you run the names command, you may note that our desired variables are in columns 2, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15.

names(WELLS_CV_NO3.sf)
WELLS_CV_NO3.sf <- dplyr::select(WELLS_CV_NO3.sf, 2, 4:5, 7:15)

Here an alternative approach if you prefer selecting variables by name instead of number:

WELLS_CV_NO3.sf <- dplyr::select(WELLS_CV_NO3.sf, c("GM_Results", "well id", "chemical id", "UNIT", "WELL_DEPTH_FT", "data source", "WELL_TYPE", "TOP_OF_SCREEN", "BOTTOM_OF_SCREEN", "DATE", "county", "geometry"))

Notice our dataset called WELLS_CV_NO3.sf in our Global Environment now has 12 variables

Which method of selecting variables do you like better? Why?

7. Rename variables

It’s useful to rename variables, especially shortening names, removing spaces, and making your variables consistent across datasets:

WELLS_CV_NO3.sf <- rename(WELLS_CV_NO3.sf, 
                              RESULTS = `GM_Results`,
                              WELL_ID = `well id`,
                              CHEMICAL_ID = `chemical id`,
                              DATA_SOURCE = `data source`,
                              COUNTY = `county`,
                              GEOMETRY = `geometry`)

Notice which side of the equation the “new” and “old” names go on. Also what do the `` characters do in our syntax?

8. Add new variables (mutate)

Let’s add two new variables. The first one adds a column that specifies the MCL for nitrate (10 mg/L), the second then calculates what the percentage of the MCL the well test result represents.

WELLS_CV_NO3.sf <- mutate(WELLS_CV_NO3.sf, MCL = as.numeric(10))
WELLS_CV_NO3.sf <- mutate(WELLS_CV_NO3.sf, PERC_MCL = as.numeric(100*RESULTS/MCL))

Let’s create an indicator variable that shows which results were above the MCL:

WELLS_CV_NO3.sf <- mutate(WELLS_CV_NO3.sf, ABOVE_MCL = ifelse(RESULTS > MCL, 1, 0))
WELLS_CV_NO3.sf$ABOVE_MCL <- as.factor(WELLS_CV_NO3.sf$ABOVE_MCL)

# Or, to be even more efficient, you can combine two lines of code into one:

WELLS_CV_NO3.sf <- mutate(WELLS_CV_NO3.sf, ABOVE_MCL = as.factor(ifelse(RESULTS > MCL, 1, 0)))

9. Manage missing values appropriately

When you work with data in R, missing values are often shown as ‘na’ in your data. For numeric data, standard values such as ‘-99’ might also be used.

Let’s start by counting the number of well tests that do not have any reported results, and dropping those that do not. Counting again should now result in zero missing RESULTS values.

sum(is.na(WELLS_CV_NO3.sf$RESULTS))
WELLS_CV_NO3.sf <- drop_na(WELLS_CV_NO3.sf, RESULTS)
sum(is.na(WELLS_CV_NO3.sf$RESULTS))

What are the benefits and drawbacks of the above approach?

Next, we will count the instances of missing TOP_OF_SCREEN and BOTTOM_OF_SCREEN values and replace them with “-99”

sum(is.na(WELLS_CV_NO3.sf$TOP_OF_SCREEN))
WELLS_CV_NO3.sf <- replace_na(WELLS_CV_NO3.sf, list(TOP_OF_SCREEN = -99))

sum(is.na(WELLS_CV_NO3.sf$BOTTOM_OF_SCREEN))
WELLS_CV_NO3.sf <- replace_na(WELLS_CV_NO3.sf, list(BOTTOM_OF_SCREEN = -99))

10. Re-write steps 4-9 using pipes and efficient coding

Check out how I used the pipe (%>%) and consolidated the code. the pipe function allows my code to gracefully transition between functions without having to have extra ouputs or naming an object each time. To attempt this you will have to reread in the sf object

WELLS_CV_NO3_v2.sf <- st_as_sf(WELLS_CV_NO3, coords = c("GM_Longitude", "GM_Latitude"), 
                             crs = 4269, agr = "constant") # Convert to sf, crs = NAD 1983
WELLS_CV_NO3_v2.sf <- st_transform(WELLS_CV_NO3_v2.sf, crs = 3310)

WELLS_CV_NO3_v2.sf <- WELLS_CV_NO3_v2.sf %>%
  distinct(.keep_all = TRUE) %>%
  filter(county  == "stanislaus") %>%
  dplyr::select(2, 4:5, 7:15) %>%
  rename(RESULTS = `GM_Results`,
    WELL_ID = `well id`,
    CHEMICAL_ID = `chemical id`,
    DATA_SOURCE = `data source`,
    COUNTY = `county`,
    GEOMETRY = `geometry`) %>%
  mutate(MCL = as.numeric(10),
         PERC_MCL = as.numeric(100*RESULTS/MCL),
         ABOVE_MCL = as.factor(ifelse(RESULTS > MCL, 1, 0)))  %>%
  drop_na(RESULTS) %>%
  replace_na(list(TOP_OF_SCREEN = -99)) %>%
  replace_na(list(BOTTOM_OF_SCREEN = -99))

Read through this code and write in simple language what each line does. What is the percentage of well test results that were above the MCL in San Joaquin, Merced, and Stanislaus Counties?

11. Organize your RScript and save

You might want to create some categories like these:

# Set up

# Data import

# Data cleaning

Applying what you’ve learned

Find a spatial or non-spatial dataset that you may use in your final project.
Assess the status of the data; how does it look in terms of Wickham’s tidy criteria?
Tidy your data, applying at least four of the the following functions, as needed: distinct, dplyr::filter, dplyr::select, rename, and mutate.
Do you have any missing values? What is your preferred strategy for these missing values?
Please upload to Camino: your final RScript and a short (1 paragraph) description of how you made your dataset more tidy.

ENVS 117 Module 3 Lab - Tidying data

Jake Dialesandro