Coming out of depression and into birding

I've heard my story from so many others in the broad grad student/academic community on Twitter. After four years in undergraduate relatively close to home, family, and friends; I moved one thousand miles west to Michigan. I had worked previously as an RA living in residence halls, and I went from that to a studio apartment alone but close to school. Lacking any connections or a support network in my new town, my mood quickly destabilized. I would be told later that while I had built a strong support network in undergraduate, I may have been pre-disposed to depression. I'm still not sure.

I can't say I was depressed before college, and so I was unprepared to recognize or deal with this change. In an effort to assimilate into a community, I joined the swing dancing community in West Michigan (an activity from my undergraduate). Something just didn't feel right though, and I didn't make connections the way I had in undergraduate. The starkest contrast between grad school and undergraduate was my daily interactions with friends, classmates, and peers contrasted with my new office: where I might see two people daily. Like anyone who has dealt with depression knows, the effort to reach out - even to join group coffee runs - began to be sapped away. As summer 2017 came and went, the winter blues had come to stay.

Perhaps one of the only things that kept me above the surface was fostering greyhounds. For a grad student, it matched with my irregular travel schedule, as each was a commitment of a few weeks. I fostered two more, to make five total, in fall 2017. This also kept me outside, exploring the parks around me.

View this post on Instagram

Thankful for all the adventures with my two lovely fosters this semester! 🐕 Helping make grad school more bearable. And so glad for all of the greyhounds full of love, naps, and silliness ❤️

A post shared by Carson🐕 and Tyler🕺 (@longodoggorescues) on Nov 23, 2017 at 1:17pm PST

I had a bird feeder on my porch, and for much of the time I spent exhausted on my couch that fall and winter, I could look outside and see woodpeckers, cardinals, and chickadees. I also had a very noisy nest of house sparrows living in the roof right next to my bedroom.

As I got more variety, like the little Junco shown above, I felt some of my curiosity perked. The Merlin App (link) helped me with my ID skills, and by that spring I was able to recognize a Scarlet Tanager when it visited me!

Something else changed too over that winter. In February 2018, I quit drinking - a habit that had become destructive and would isolate me to the couch or my floor, alone - and after a year of in-person therapy I sought medical help from the University psychiatry office.

It took a few months to get the dosing right, but in April my mood - plagued before by dramatic swings (yes, I kept the data!) - dramatically improved and stabilized. Thankfully, this improvement came a week before my MS defense. I still credit my seeking help with giving me enough willpower to write/finish my MS thesis in the first place. 

Like a case of bad congestion, after I started seeing improvements suddenly the woods would be full of light and flitting birds and songs. I moved to Waterloo to begin my PhD. I'm still pursuing the same course of study, in the earth and natural sciences. Fortunately, my work offers the chances to travel and see little details in rocks, trees, and landscapes. I suppose birding fits naturally into this way of seeing the world. 

I've since adopted a rescue greyhound, who is very tolerant of my habit of stopping on walks to just listen. Begrudgingly, we're outside all year round, but bundled up. I'm taking Vitamin D too in the winter, but birding gives one more interesting reason to enjoy (if only a little) the outdoors in freezing temperatures. 

It's just been a year since starting medication and quitting drinking. I'm again in a new place, and new communities, but it's so much better than last time. Of course, this speaks volumes of my new lab group, and academic adviser, and my canine companion. I find myself just as passionate about my research, and I've made tremendous progress (so I'm told) in just two terms. But the birding itch now has to be scratched. As the weather warms and sun sets later, I'm outside for a good portion of the mornings and evenings, and my thoughts are free to roam between my research, what flowers to plant this spring, and naming the bird songs that ring through the neighborhood woodlot. I can't help but find it ironic that even my new hobby is all about science, as contributions of photos and observations to eBird represent one of the world's largest citizen science projects.

I look forward to sharing more parts of my story. Happy birding! 🐦😊

Population of Canada within 100 miles of the US border, with Open Data and R!

Introduction

This tutorial uses the sf package to perform geoprocessing, and ggplot to visualize the data. Additionally, we practice downloading open-source data from the internet.

Adjust Settings, Load Packages

knitr::opts_chunk$set(echo = TRUE)
for(p in c("ggplot2","sf","dplyr","plyr","raster")){library(p,character.only = T)}
directory=".../Layers_Canada/CanadaPopulationDensity/"
theme_set(theme_bw()+
            theme(
              panel.background = element_blank(),axis.ticks = element_blank(),
              legend.background = element_blank(),panel.grid = element_line(colour = 0)
            ))

Download Open-Source data from open.canada.ca

setwd(directory)
if(!("CAdata.zip" %in% list.files())){
  #https://open.canada.ca/data/en/dataset/5be03a46-8504-40a7-a96c-af195bae0428
  download.file(url = "http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/gecu000e11a_e.zip",
                dest=paste0(directory,"CAdata.zip"))
}
if(!("CAcensusdata" %in% list.files())){unzip(zipfile = "CAdata.zip",exdir = "CAcensusdata")}
#https://open.canada.ca/data/en/dataset/ece81c43-aa4e-41ef-86c2-3835eb5aa95c
download.file(url = "http://www12.statcan.gc.ca/census-recensement/2016/dp-pd/hlt-fst/pd-pl/Tables/CompFile.cfm?Lang=Eng&T=701&OFT=FULLCSV",dest=paste0(directory,"CAcensusdata/CensusData.csv"))
if(!("Borders.zip" %in% list.files())){
  #http://thematicmapping.org/downloads/world_borders.php
  download.file(url = "http://thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip",dest=paste0(directory,"Borders.zip"))
}
setwd(setwd(directory))
if(!("WorldBorders" %in% list.files())){unzip(zipfile = "Borders.zip",exdir = "WorldBorders")}

Inspect the Canadian Census data from 2016.

We see that population within census divisions is log-normally distributed. We can also confirm that the population density calculations within each census division match the population and area values. The population of Canada in 2016 was 35.2 million people.

census=read.csv(paste0(directory,"CAcensusdata/CensusData.csv"))
ggplot(data=census)+
  geom_histogram(aes(x=Population..2016))+
  scale_x_log10()

# Confirm that Population Density calculations match areas
ggplot(data=census)+
  geom_point(aes(y=Population..2016,x=Land.area.in.square.kilometres..2016*Population.density.per.square.kilometre..2016))+
  geom_abline(intercept=0,slope=1)

print(paste(round(sum(census$Population..2016,na.rm=TRUE)*10^-6,1),"million people"))

## [1] "35.2 million people"

Inspect GIS Data

proj="+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

CAshp=st_read(paste0(directory,"CAcensusdata/gcd_000e11a_e.shp"))

## Reading layer `gcd_000e11a_e' from data source `...\Layers_Canada\CanadaPopulationDensity\CAcensusdata\gcd_000e11a_e.shp' using driver `ESRI Shapefile'
## Simple feature collection with 293 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -141.0181 ymin: 41.71096 xmax: -52.61941 ymax: 83.1355
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs

CAshp=st_transform(CAshp,proj)
CAshp=dplyr::left_join(CAshp,census,by=c("CDUID"="Geographic.code"))

ggplot()+
  geom_sf(data=CAshp)+
  ggtitle("Canada Census Divisions")

borders=st_read(paste0(directory,"WorldBorders/TM_WORLD_BORDERS-0.3.shp"))

## Reading layer `TM_WORLD_BORDERS-0.3' from data source `...\Layers_Canada\CanadaPopulationDensity\WorldBorders\TM_WORLD_BORDERS-0.3.shp' using driver `ESRI Shapefile'
## Simple feature collection with 246 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.6236
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

USA=borders[borders$NAME=="United States",c("NAME","geometry")]
USA=st_transform(USA,proj)
st_crs(USA)$units

## [1] "m"

USA_buf=st_buffer(USA,dist=160934) # 100 miles = 160934 meters

ggplot()+
  geom_sf(data=USA,fill="gray",col="transparent")+
  geom_sf(data=USA_buf,fill="transparent",col=2)+
  ggtitle("US Border and 100 mile buffer")

Inspect Population Density Data

We can see that the region within 100 miles of the US border is highly populated, but there are cities like Calgary, Edmonton, and Saskatoon that are north of this boundary.

ggplot()+
  geom_sf(data=CAshp,aes(fill=log10(Population.density.per.square.kilometre..2016)),col="transparent")+
  geom_sf(data=USA_buf,col=2,fill="transparent")+
  coord_sf(xlim = extent(CAshp)[1:2],ylim = extent(CAshp)[3:4])+
  ggtitle("Canadian Census Divisions and Population Density")+
  scale_fill_gradient(name="Pop_dens",low="darkblue",high="red",na.value="darkblue",
                      limits=c(-1.2,3.5),breaks=c(-1:3),labels=(10^c(-1:3)))

Use st_intersection to extract census tracts

We use st_intersection to intersect the buffered US borer polygon and Canadian census tracts. We assume population is evenly distributed within census divisions, allowing us to adjust the size of intersected tracts for their new clipped size. Inspecting their new calculated areas, we see that many tracts are completely within the 100 mile buffer, but others have a much lower area.
Our final calculation finds that 27.3 million people live within 100 miles of the US border, or 77.7% of the population.

intersect=st_intersection(CAshp,USA_buf)
plot=ggplot()+
  geom_sf(data=CAshp,fill="gray",col="transparent")+
  geom_sf(data=intersect,aes(fill=log10(Population.density.per.square.kilometre..2016)),col="transparent")+
  ggtitle("Population of Canada within 100 miles of the US border")+
  scale_fill_gradient(name=expression(paste("Persons km",NA^-2)),
                      low="darkblue",high="red",na.value="darkblue",
                      limits=c(-1.2,3.5),breaks=c(-1:3),labels=(10^c(-1:3)))+
  theme(legend.position = c(0.85,0.8),axis.ticks = element_blank(),
        axis.text = element_blank(),axis.title = element_blank())
ggsave(plot,filename = paste0(directory,"BufferPopDens.png"),width = 5,height=4,dpi=300,units="in")
plot

intersect=intersect%>%mutate(area_m2 = st_area(.) %>% as.numeric())
intersect$area_km2=intersect$area_m2*(1000^-2)

ggplot(data=intersect)+
  geom_point(aes(y=area_km2,x=Land.area.in.square.kilometres..2016))+
  scale_y_log10()+scale_x_log10()+
  geom_abline(intercept=0,slope=1)

intersect$pop_2016_intersect=intersect$area_km2*intersect$Population.density.per.square.kilometre..2016
paste0(round(sum(intersect$pop_2016_intersect)*10^-6,1)," million people")

## [1] "27.3 million people"

print(paste0(round(100*sum(intersect$pop_2016_intersect,na.rm=TRUE)/sum(census$Population..2016,na.rm=TRUE),1),"% of population within 100 miles of the US border"))

## [1] "77.7% of population within 100 miles of the US border"

What about Canada's population in relation to some US city? Seattle?

lon=seq(-140,-45,2)
df <- data.frame(lon=c(lon,rev(lon)),lat=c(rep(47.6,length(lon)),rep(30,length(lon))))
df <- rbind(df, df[1,])
seattle <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))), crs = 4326)
seattle=st_transform(seattle,proj)

intersect=st_intersection(seattle,CAshp)
intersect=intersect%>%mutate(area_m2 = st_area(.) %>% as.numeric())
intersect$area_km2=intersect$area_m2*(1000^-2)
intersect$pop_2016_intersect=intersect$area_km2*intersect$Population.density.per.square.kilometre..2016


ggplot()+
  geom_sf(data=CAshp,fill="gray",col="transparent")+
  geom_sf(data=intersect,aes(fill=log10(Population.density.per.square.kilometre..2016)),col="transparent")+
  geom_sf(data=seattle,col=2,fill="transparent")+
  coord_sf(xlim = extent(CAshp)[1:2],ylim = extent(CAshp)[3:4])+
  ggtitle("Population of Canada south of Seattle, Washington")+
  scale_fill_gradient(name=expression(paste("Persons km",NA^-2)),
                      low="darkblue",high="red",na.value="darkblue",
                      limits=c(-1.2,3.5),breaks=c(-1:3),labels=(10^c(-1:3)))+
  theme(legend.position = c(0.85,0.8),axis.ticks = element_blank(),
        axis.text = element_blank(),axis.title = element_blank())

paste0(round(sum(intersect$pop_2016_intersect)*10^-6,1)," million people")

## [1] "23.3 million people"

print(paste0(round(100*sum(intersect$pop_2016_intersect,na.rm=TRUE)/sum(census$Population..2016,na.rm=TRUE),1),"% of population south of Seattle, Washington State"))

## [1] "66.2% of population south of Seattle, Washington State"