<jupyter_start><jupyter_text>Session Prep 1. Run the cell below to install two libraries we'll be using to look at geospatial data: GeoPandas and OpticalRS.2. Visit the [Github repo](https://github.com/rrcarlson/Intro_python_geo) and download/clone the files in the Data folder (shapefiles and raster).<jupyter_code>#@title
# Install OpticalRS. This also installs GeoPandas because GeoPandas is a dependency of OpticalRS
!apt-get install software-properties-common python-software-properties > /dev/null
!add-apt-repository ppa:ubuntugis/ppa -y > /dev/null
!apt-get update > /dev/null
!apt-get install -y --fix-missing python-gdal gdal-bin libgdal-dev > /dev/null
!pip2 install OpticalRS > /dev/null<jupyter_output><empty_output><jupyter_text>Intro to Python ![Python image](https://github.com/rrcarlson/Intro_python_geo/blob/master/images/python.jpg?raw=1)Python is a rich, versatile language with lots of prolific contributors. So where to start? For reference, below are some Python tutorials that cover a range of Python libraries. We won't spend too much time today on basic syntax, etc., since that's well covered elsewhere and kinda boring. As a Python newbie, I've found that using Pandas/GeoPandas and NumPy have introduced me organically to other libraries, and have allowed me to see results quickly for basic geospatial tasks--so we'll start there!* [Python basics](https://github.com/sinkovit/PythonSeries/blob/master/Python%20basics.ipynb): How to create lists, dictionaries, iterators, etc.* [Python Markdown and LaTex](https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet): Documentation and formatting* Matplotlib: Data visualization* [Pandas](https://pandas.pydata.org/pandas-docs/stable/10min.html): Data wrangling and preparation using two types of data objects: 1) DataFrames, 2) Series.* [GeoPandas](https://automating-gis-processes.github.io/2016/Lesson2-overview-pandas-geopandas.html): Basically the geospatial version of Pandas. Adds geospatial functionality to the DataFrame (Geo + DataFrame = GeoDataFrame).* [NumPy](https://www.datacamp.com/community/tutorials/python-numpy-tutorial): Multi-dimensional data arrays for manipulating large datasets. There are some functions for linear algebra and statistics within NumPy but it's also used by SciPy, matplotlib, and pandas for scientific computing.* [SciPy](http://www.randalolson.com/2012/08/06/statistical-analysis-made-easy-in-python/): Scientific computing and statistical analysis. Uses NumPy, Matplotlib, and Pandas, and is common in earth sciences, astronomy, and oceanography.* [Shapely](http://toblerity.org/shapely/manual.html): Spatial data model for planar features (points, curves, and surfaces).* Skikit-learn: regression, clustering, and classification algorithms* Various [cool tricks](https://community.modeanalytics.com/python/)  Conservation in the Dominican Republic ![Reef image](https://github.com/rrcarlson/Intro_python_geo/blob/master/images/01_coral_reef_bright_spot.jpg?raw=1)In this tutorial, we'll be testing out two Python libraries, **GeoPandas** and **NumPy**, that are useful in GIS. Our objectives: to filter and analyze vector and raster data for conservation areas in the Dominican Republic.We'll use two datasets:* **Conservation Areas of the Dominican Republic**: Vector data delineating conservation area boundaries, from the [World Database of Protected Areas (WDPA)](https://www.iucn.org/theme/protected-areas/our-work/world-database-protected-areas).* **Gridded Bathymetric Data**: Raster data (30 arc-second resolution) on bathymetry from the [British Oceanographic Data Center](https://www.bodc.ac.uk/).Our goals are:* **Filter the WDPA database** to find only conservation areas for the marine environment (Marine Protected Areas)* **Clip rasters** to MPA boundaries* **Analyze raster data** (bathymetry) within Marine Protected Area boundaries to find    * Mean depth per MPA    * Standard deviation of depth per MPA    * Range of depth per MPA  Step 1: Load your libraries What modules do we need?* Vector data comes in a tabular format, with one geometric feature per row. So we need a library to organize geographic data in a DataFrame.**That's Pandas** (more specifically, **GeoPandas**).* Raster data is a little trickier. It comes to us in a grid of pixel values, sometimes in multiple bands (3+ dimensions). We need a library to organize arrays/matrices of values in *n* dimensions. **That's NumPy**.* Numpy isn't explicitly designed to handle raster data. We need the help of another library to convert our raster files into array format. **That's [OpticalRS](https://github.com/jkibele/OpticalRS)**, created by NCEAS' very own <NAME>! OpticalRS borrows a bit of code from [an old version of rasterstats](https://github.com/perrygeo/python-rasterstats/releases/tag/0.5) to do the actual raster subsetting.* We'll want to plot stuff, so we need a library for data visualization. That's **matplotlib**.![arrays image](https://github.com/rrcarlson/Intro_python_geo/blob/master/images/arrays.png?raw=1)<jupyter_code>import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib
from OpticalRS import *<jupyter_output><empty_output><jupyter_text>Step 2: Prep your vector data  Step 2.1 Read in your vector data: GeoPandas Reading in tabular data is pretty simple. Just remember:* DataFrame = Pandas = `pd.read_csv()`* GeoDataFrame = DataFrame with geometry attribute = GeoPandas = `gpd.read_file()`<jupyter_code># Loading my local files into Colaboratory
from google.colab import files
uploaded = files.upload()
# Define your vector filepath
conserve_fp = "WDPA_polygons.shp"

# Read vector data in as GeoDataFrame
conserve = gpd.read_file(conserve_fp)<jupyter_output><empty_output><jupyter_text>Our first preprocessing step is to check on the data's projection. Projections can be tricky and I'm not going to delve too far into this topic here, but below are some basic steps for identifying and changing your data's projection (more details can be found [here](http://geopandas.org/projections.html)).<jupyter_code># Check the Coordinate Reference System (CRS)
conserve.crs
# Change your CRS, for fun. Then change it back (because we like EPSG = 4326).
conserve = conserve.to_crs(epsg=3857)
conserve = conserve.to_crs(epsg=4326)<jupyter_output><empty_output><jupyter_text>Exercise 2.1<jupyter_code># Read in your WDPA shapefile as a GeoDataFrame.

# Check the projection of your data.<jupyter_output><empty_output><jupyter_text>Step 2.2 Explore your vector data: GeoPandas We can use some basic exploratory functions in GeoPandas to query our data.<jupyter_code># Look at the first few features in our attribute table
conserve.head()
# Look at all of the attributes in your dataset
conserve.info()<jupyter_output><empty_output><jupyter_text>There are 143 data objects with 29 attributes, of various data types (object, float, integer). In the pandas world, data types break down like this:* dtype('O') = object = can include string and other geopandas object types* dtype('int64') = integer* dtype('float64') = float<jupyter_code># What is the total conserved area in the DR? Conserved area = "GIS_AREA" (we'll discuss how to use the inherent geometry of the GeoDataFrame to do this later).
area_sum = conserve.GIS_AREA.sum()
area_sum
# What is the largest area conserved?
area_max = conserve.GIS_AREA.max()
area_max
# What is the smallest area conserved?
area_min = conserve.GIS_AREA.min()
area_min
# What are the unique values for the "MARINE" attribute?
conserve['MARINE'].unique()
# How many different areas are "Marine" (2), "Terrestrial" (0), or both (1)?
type_park = conserve['MARINE'].value_counts()
type_park
# Write all of this info into a string
"The total conserved area is {} square meters and there are {} MPAs, {} terrestrial areas, and {} mixed areas".format(area_sum, type_park["2"], type_park["0"], type_park["1"])<jupyter_output><empty_output><jupyter_text>Exercise 2.2<jupyter_code># What are the options for park designation in English ("DESIG_ENG")?

# How many different areas are designated as "Wildlife Refuge"?<jupyter_output><empty_output><jupyter_text>For a complete list of methods available for exploring GeoDataFrames, see [here](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html).  Step 2.3 Plot your vector data: GeoPandas GeoPandas offers some basic plotting methods for quick visualization. You can create more advanced plots (with customizable layering, formatting) using the Matplotlib library. Matplotlib can be a bit confusing at first, so here's a [great tutorial](https://github.com/matplotlib/AnatomyOfMatplotlib/blob/master/AnatomyOfMatplotlib-Part1-Figures_Subplots_and_layouts.ipynb) if you want to learn more. As in R, you can also create interactive maps through [Leaflet](https://automating-gis-processes.github.io/2016/Lesson5-interactive-map-folium.html) and [Bokeh](https://automating-gis-processes.github.io/2016/Lesson5-interactive-map-bokeh.html).<jupyter_code># Plot conservation areas
conserve.plot()
# Plot conservation areas as choropleth by park type
conserve.plot(column = "DESIG_ENG", cmap = "Paired")
# Make it red! Color palettes: https://matplotlib.org/users/colormaps.html
conserve.plot(column = "DESIG_ENG", cmap = "hot")
# Plot conservation areas where the area has a Use Management Plan.
conserve[conserve['MANG_PLAN'] != "No"].plot()<jupyter_output><empty_output><jupyter_text>Exercise 2.3<jupyter_code># Plot conservation areas as a choropleth by the "MARINE" attribute
conserve.plot(column = "MARINE")
# Make it blue (obviously)! https://matplotlib.org/users/colormaps.html
conserve.plot(column = 'MARINE', cmap="Blues")<jupyter_output><empty_output><jupyter_text>Step 2.4 Filter your vector data: GeoPandas **Looking above at our "To Do" list, our first task was this:** "Filter the WDPA database to find only conservation areas for the marine environment (Marine Protected Areas)."In R-speak, we're using Python to `dplyr::filter`. To subset data in Python, we use the `loc` and `iloc` indexing tools. Documentation can be found [here](http://pandas.pydata.org/pandas-docs/version/0.17/indexing.html) under "Indexing and Selecting Data" (basically the Python version of R's Data Wrangling cheat sheet).The basics are this:* `iloc` selects data on a known index (integer value). Example: `iris.iloc[7]` would select the 8th row of data.* `loc` selects data on a known label. Example `iris.loc['sepal_length'>7]` selects data with sepal length over 7 centimeters.Let's try it out:<jupyter_code># Take a look at our first few lines of data
conserve.head(3)
# Select rows 0 through 2, column 4
conserve.iloc[0:3, 3]
# Select index value 3 THROUGH 6 for name and metadata ID
conserve.loc[3:6, ['NAME','METADATAID']]
# Select indices 3 AND 6 for name and metadata ID
conserve.loc[[3,6], ['NAME','METADATAID']]
# Filter DataFrame for only those areas with the name "Jaragua"
conserve.loc[conserve['NAME'] == "Jaragua"]<jupyter_output><empty_output><jupyter_text>The WDPA Data Dictionary that Marine Protected Areas are coded by the ordinal `MARINE` attribute, defined as follows:* 0 = 100% terrestrial PA* 1 = Coastal: marine and terrestrial PA* 2 = 100 % marine PALet's subset our data by *only* those rows that represent 100% MPAs.<jupyter_code># Select only 100% marine protected areas
mpa = conserve.loc[conserve['MARINE'] == '2']
len(mpa)
mpa[['MARINE', 'NAME']]<jupyter_output><empty_output><jupyter_text>Our resulting DataFrame contains vector data for 7 MPAs in the Dominican Republic.  Before we move on... **Can you "pipe" in Python?** Yes and no. You won't see any "%>%" business in Python, but Python syntax allows you to string together multiple operations. The Python mantra is that "everything is an object", meaning everything (strings, dataframes, lists, functions, even modules) can be assigned to a variable or passed as an argument to a function.Confusing? Take a look at our code above. We're running multiple operations through our "conserve" geodataframe, like: 1) reprojecting it, 2) filtering it, and 3) plotting it. In R, "piping" might look something like: `conserve %>%`> `st_transform(4326) %>%`    > `dplyr::filter(MARINE == "2") %>%`   > `dplyr::select(MARINE, NAME, geometry) %>%`> `plot()`    In Python, "piping" works out like this:`conserve.to_crs(epsg=4326).loc[conserve['MARINE']=='2'][['MARINE','NAME','geometry']].plot()`Since "everything is an object" in Python, the result of each operation can be worked on by subsequent methods. You can see that this line of code simply combines the functions we used in Step 2.2. Let's see if it works.<jupyter_code>conserve.to_crs(epsg=4326).loc[conserve['MARINE'] == '2'][['MARINE','GIS_AREA','geometry']].plot()<jupyter_output><empty_output><jupyter_text>Woohoo!  Also before we move on... **You can also perform call kinds of operations using the 'geometry' column (every GeoDataFrame has one).**Above, we used the dataset's "GIS_AREA" attribute to run summary statistics. This is because the WDPA dataset uses its own methodology and decision factors for calculating polygon areas. Most datasets won't be this finicky.The "geometry" column contains a Shapely geometry object (point, polygon, multi-point, multi-polygon, etc.) and is usually automatically detected when GeoPandas reads in your shapefile. But be careful! The geometry column won't always be called "geometry". Use the "name" method to double check. If GeoPandas chooses the wrong column for geometry, use the `set_geometry` method to specify which column to use.<jupyter_code># Call your geometry attribute and find its Shapely geometry type
conserve.geometry.head()
# Check which column GeoPandas has chosen as the "geometry" attribute
conserve.geometry.name
# Change if necessary (not necessary in our case, so this code is redundant)
conserve.set_geometry("geometry")
conserve.geometry.name<jupyter_output><empty_output><jupyter_text>Step 3: Prep your raster data  Step 3.1 Read in your raster data: Numpy and OpticalRS **Let's turn our attention to our raster data.** One of my favorite parts of Python is its ability to manipulate raster data. By turning a raster into an object called a NumPy array, you can directly access the matrix of values within raster pixels, analyzing and filtering raster data in [all kinds of ways](https://docs.scipy.org/doc/numpy/user/quickstart.html) within the expansive SciPy library. Arrays are also compatible with GeoPandas, so you don't have to worry about converting your vector data into SpatialPolygonsDataFrames (or whatever) to play well with rasters.To turn our raster file into a numpy array, we'll use a library called OpticalRS by <NAME> found [here](https://github.com/jkibele/OpticalRS).<jupyter_code>%pylab inline
from OpticalRS import *<jupyter_output><empty_output><jupyter_text>Our first step is to read our raster file into OpticalRS. The `RasterDS()` function reads in a file and converts it into a raster-specific object called a 'Raster Dataset'. The `band_array` attribute of the 'Raster Dataset' is a NumPy array of your data.<jupyter_code># Loading my local files into Colaboratory
from google.colab import files
uploaded = files.upload()
# Define your GeoTiff filepath
bathy_dr_fp = "bathy_caribb_4326.tif"
        
# Turn your GeoTiff into a "Raster Dataset" object.
rds = RasterDS(bathy_dr_fp)

# Convert the "Raster Dataset" object into a numpy array
dominican_arr = rds.band_array<jupyter_output><empty_output><jupyter_text>If you have lots of rasters to read in at once, define your root directory and run a loop through it```in_data_dir = "/home/carlson/tutorial_files/"filenames = []for p in os.listdir(in_data_dir):    if p.endswith('.tif'):        fp = os.path.join(in_data_dir, p)        filenames.append(fp)``` Let's check out the properties of our array and make sure the conversion worked<jupyter_code># Check array shape
dominican_arr.shape
# Find mean depth (pixel value) in your raster
dominican_arr.mean()
# Check size of array
dominican_arr.nbytes
# Create a quick plot
imshow(dominican_arr.squeeze(), cmap ='viridis')<jupyter_output><empty_output><jupyter_text>Looks good!  Exercise 3.1<jupyter_code># Follow the steps to load in your raster data as a NumPy array.
# Find the max bathymetry value in your array (use NumPy documentation: https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.statistics.html)<jupyter_output><empty_output><jupyter_text>Step 4: Analyze bathymetry by MPA ("Zonal Statistics")  Step 4.1 Apply bathymetry to one "zone" Running statistics on our array/raster is easy, but not terribly meaningful because my bounding box is pretty arbitrary. We need to clip our raster to MPAs in order to assess the bathymetry values within our "zones" of interest.  This is called "zonal statistics" (aka "cookie cutter statistics"). ![Dwayne the rock image](https://github.com/rrcarlson/Intro_python_geo/blob/master/images/the_rock.png?raw=1)Let's start by clipping our "Raster Dataset" object to MPA geometries. This is easy in OpticalRS using the `geometry_subset()` function. `geometry_subset` takes a geometry object (like the 'geometry' attribute of a GeoDataFrame row) and uses it as a bounding box (box-ish) for our "Raster Dataset", returning a NumPy array.Since `geometry_subset` takes just one geometry, we'll1. Dissolve all parks designated as "Marine" into one geometry.2. Use this geometry as the "cookie cutter" for our raster.<jupyter_code># Create one geometry for ALL marine protected areas (Remember: MPAs are coded as '1')
mpas_dissolved = conserve.dissolve(by='MARINE').loc['2']

# Clip your array to the shape of MPAs (dissolved)
mpa_array = rds.geometry_subset(mpas_dissolved.geometry)

# Plot the clipped array
imshow(mpa_array.squeeze())<jupyter_output><empty_output><jupyter_text>Great! Now we have an array of bathymetry values clipped to MPA boundaries. We can use [any statistical methods](https://docs.scipy.org/doc/numpy-1.14.0/reference/index.html) available within NumPy on this array.<jupyter_code># Find the mean, min, max, and range of bathymetry values for MPAs.
[mpa_array.mean(), mpa_array.min(), mpa_array.max(), mpa_array.ptp()]<jupyter_output><empty_output><jupyter_text>Let's tighten this up by placing these commands in this same function. **Note:** There are lots of preexisting functions for doing "Raster Stats" in Python and R, but they're prescriptive in the types of statistics they return. By converting your raster into an array, you can apply anything in the **NumPy** universe to your data!<jupyter_code># Define a function encompassing all of the commands above.
# Clip Raster Dataset to vector object g.
# Find mean, min, max, and range values of clipped numpy array.
# Fetch the MPA ID from "g" and place alongside summary statistics.
def stats_bathy(g, rds):
    mpa_arr = rds.geometry_subset(g.geometry)
    mean_bathy = mpa_arr.mean()
    min_bathy = mpa_arr.min()
    max_bathy = mpa_arr.max()
    range_bathy = mpa_arr.ptp()
    mpa_id = g['WDPA_PID']
    return mpa_id, mean_bathy, min_bathy, max_bathy, range_bathy
# Try the function on g = mpas_dissolved
stats_bathy(mpas_dissolved, rds)<jupyter_output><empty_output><jupyter_text>Exercise 4.1<jupyter_code># Follow the steps above to clip your raster data to MPA boundaries.
# Define a new function calculating the standard deviation and coefficient of variation of MPA geometries<jupyter_output><empty_output><jupyter_text>Step 4.2: Apply bathymetry to many "zones" Now that our personal "Raster Stats" are bundled neatly in one function, we can iterate through multiple geometries, for example, each individual MPA (undissolved). The easiest way to do this is to use `apply`, which operates like `apply` in R. Our steps are:* Create an anonymous function using `lambda` and run each geometry (row) in our DataFrame through that function. * `apply` basically **points each geometry in our DataFrame**  to fill `g` in `stats_bathy(g, rds)`<jupyter_code># Define a new function that uses 'lambda' to iterate stats_bathy through multiple DataFrame rows (MPA geometries).
def stats_bathy_multiple(dataframe, rds):
    series_mpas = dataframe.apply(lambda g: stats_bathy(g, rds), axis=1)
    return series_mpas
# Try out the function on our original 'mpa' GeoDataFrame
series_mpas = stats_bathy_multiple(mpa, rds)
series_mpas<jupyter_output><empty_output><jupyter_text>We now have summary statistics for all MPAs in our `conserve` dataset!  Extra steps The data above isn't in the most usable format. Our function returns a series of tuples, with each tuple representing a different MPA geometry. If I want to get all of this back into a DataFrame, I can split each tuple into separate values (attributes). There are several ways to do this, but one easy method is `vstack()`, which splits *thing* into multiple *sub-things* and stacks them vertically to make an array. I can then insert this data into a DataFrame.<jupyter_code># Convert your series of tuples into a DataFrame using 'vstack'
dataframe_mpas = pd.DataFrame(np.vstack(series_mpas))
dataframe_mpas
# Rename the columns of your DataFrame
dataframe_mpas = dataframe_mpas.rename(columns = {0:'WDPA_PID', 1:'mean_bathy', 2: 'min_bathy', 3: 'max_bathy', 4: 'range_bathy'})
dataframe_mpas<jupyter_output><empty_output><jupyter_text>I can then save this data back out to a CSV. Saving tabular data is as simple as:* DataFrame = Pandas = `pd.read_csv()` = `pd.to_csv()`* GeoDataFrame = DataFrame with geometry attribute = GeoPandas = `gpd.read_file()` = `gpd.to_file()`<jupyter_code># Save DataFrame to CSV
dataframe_mpas.to_csv("python_ecodatascience.csv")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09360424844146849, 'fraction_numerical': 0.008312168090510275, 'mean_word_length': 5.6206053194741665, 'pattern_counts': {'":': 0, '<': 103, '<?xml version=': 0, '>': 118, 'https://': 25, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '214852', 'license': 'Apache-2.0', 'n_tokens_mistral': 6855, 'n_tokens_neox': 6399, 'n_words': 2514, 'path': 'copy notebooks/Copy of Intro_Python_geospatial(1).ipynb', 'repo_name': 'bsivavenu/Google_Colab_Notebooks', 'size': '21667'}