Calling R libraries from Python

In this example we will explore the Coral Reef Evaluation and Monitoring Project (CREMP) data available in the Gulf of Mexico Coastal Ocean Observing System (GCOOS) ERDDAP server.

To access the server we will use the rerddap library and export the data to Python for easier plotting.

The first step is to load the rpy2 extension that will allow us to use the R libraries.

%load_ext rpy2.ipython

The first line below has a %%R to make it an R cell. The code below specify the GCOOS server and fetches the data information for the cremp_fk_v2_1996 dataset.

For more information on rerddap please see https://rmendels.github.io/Using_rerddap.nb.html.

%%R

library('rerddap')

url <- 'https://gcoos4.tamu.edu/erddap'
data_info <- rerddap::info('fk_CREMP_Rev_DATA_v2_3_1_1996', url=url)

data_info
<ERDDAP info> fk_CREMP_Rev_DATA_v2_3_1_1996 
 Base URL: https://gcoos4.tamu.edu/erddap 
 Dataset Type: tabledap 
 Variables:  
     acceptedNameAuthorship: 
     acceptedNameUsage: 
     acceptedNameUsageID: 
     averageNumberOfPoints: 
         Range: 1045, 2625 
         Units: count 
     basisOfRecord: 
     bottomType: 
     class: 
     country: 
     datasetID: 
     datasetName: 
     depth: 
         Range: 3.0, 54.0 
         Units: m 
     eventDate: 
         Range: 1996, 1996 
     eventID: 
     family: 
     firstYear: 
         Range: 1996, 1996 
     genus: 
     geodeticDatum: 
     habitat: 
     habitatID: 
     kingdom: 
     language: 
     lastYear: 
         Range: 2000, 2015 
     latitude: 
         Range: 24.4517, 25.2953 
         Units: degrees_north 
     license: 
     locality: 
     longitude: 
         Range: -81.9195, -80.2087 
         Units: degrees_east 
     maximumDepthInMeters: 
         Range: 3.0, 54.0 
         Units: m 
     minimumDepthInMeters: 
         Range: 3.0, 54.0 
         Units: m 
     modified: 
     occurrenceID: 
     occurrenceStatus: 
     order: 
     organismQuantity: 
         Range: 0.0, 0.919800885 
         Units: percent 
     organismQuantityType: 
     ownerInstitutionCode: 
     parentEventID: 
     phylum: 
     recordedBy: 
     Samples: 
         Range: 0, 11839 
     samplingEffort: 
     samplingProtocol: 
     scientificName: 
     scientificNameAuthorship: 
     scientificNameID: 
     siteCode: 
     siteID: 
         Range: 10, 81 
     siteName: 
     specificEpithet: 
     stateProvince: 
     stationNumber: 
         Range: 101, 754 
     subRegionID: 
     taxonomicStatus: 
     taxonRank: 
     time: 
         Range: 8.283168E8, 8.283168E8 
         Units: seconds since 1970-01-01T00:00:00Z 
     transectLengthInMeters: 
         Range: 19.0, 26.0 
         Units: m 
     type: 
     vernacularName: 
     waterBody: 

By inspecting the information above we can find the variables available in the dateset and use the tabledap function to download them.

Note that the %%R -o rdf will export the rdf variable back to the Python workspace.

%%R -o df

fields <- c(
    'depth',
    'longitude',
    'latitude',
    'organismQuantity',
    'genus',
    'habitat'
)

df <- tabledap(
    data_info,
    fields=fields,
    url=url
)
R[write to console]: info() output passed to x; setting base url to: https://gcoos4.tamu.edu/erddap

Now we need to export the R DataFrame to a pandas objects and ensure that all numeric types are numbers and not strings.

import pandas as pd


cols = ["longitude", "latitude", "depth", "organismQuantity"]
df[cols] = df[cols].apply(pd.to_numeric)

df.head()
depth longitude latitude organismQuantity genus habitat
2 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
3 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
4 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
5 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
6 9.0 -80.3782 25.1201 0.0 Acropora Hard Bottom

We can navigate to ERDDAP’s info page to find the variables description. Let’s check what is organismQuantity:

The is value of the derived information product, such as the numerical value for biomass. This term does not include units. Mean number of observed fish per species for 5 Minutes

We can see that organismQuantity has a lot of zero values, let’s remove that first to plot the data positions only where something was found.

# Filter invalid values (-999).

cremp_1996 = df.loc[df["organismQuantity"] >= 0]
cremp_1996.head()
depth longitude latitude organismQuantity genus habitat
2 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
3 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
4 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
5 6.0 -80.3475 25.1736 0.0 Acropora Hard Bottom
6 9.0 -80.3782 25.1201 0.0 Acropora Hard Bottom

What is the most common genus of Coral observed?

avg = cremp_1996.groupby("genus").mean()

avg
depth longitude latitude organismQuantity
genus
23.2375 -81.046667 24.768958 0.034680
Acropora 23.2375 -81.046667 24.768957 0.005494
Agaricia 23.2375 -81.046667 24.768957 0.000003
Cladocora 23.2375 -81.046667 24.768957 0.000000
Colpophyllia 23.2375 -81.046667 24.768957 0.004462
Dendrogyra 23.2375 -81.046667 24.768957 0.000977
Diadema 23.2375 -81.046667 24.768957 0.000007
Dichocoenia 23.2375 -81.046667 24.768957 0.000509
Dictyota 23.2375 -81.046667 24.768957 0.000000
Diploria 23.2375 -81.046667 24.768957 0.000999
Eusmilia 23.2375 -81.046667 24.768957 0.000091
Favia 23.2375 -81.046667 24.768957 0.000045
Gorgonia 23.2375 -81.046667 24.768957 0.000000
Halimeda 23.2375 -81.046667 24.768957 0.000000
Helioseris 23.2375 -81.046667 24.768957 0.000004
Isophyllia 23.2375 -81.046667 24.768957 0.000006
Lobophora 23.2375 -81.046667 24.768957 0.000000
Lyngbia 23.2375 -81.046667 24.768957 0.000000
Madracis 23.2375 -81.046667 24.768957 0.000040
Manicina 23.2375 -81.046667 24.768957 0.000007
Meandrina 23.2375 -81.046667 24.768957 0.000363
Millepora 23.2375 -81.046667 24.768957 0.006159
Montastraea 23.2375 -81.046667 24.768957 0.013597
Mussa 23.2375 -81.046667 24.768957 0.000061
Mycetophyllia 23.2375 -81.046667 24.768957 0.000144
Oculina 23.2375 -81.046667 24.768957 0.000046
Orbicella 23.2375 -81.046667 24.768957 0.033860
Palythoa 23.2375 -81.046667 24.768957 0.000000
Phyllangia 23.2375 -81.046667 24.768957 0.000000
Porites 23.2375 -81.046667 24.768957 0.003372
Pseudodiploria 23.2375 -81.046667 24.768957 0.001024
Scolymia 23.2375 -81.046667 24.768957 0.000009
Siderastrea 23.2375 -81.046667 24.768957 0.004603
Solenastrea 23.2375 -81.046667 24.768957 0.000048
Stephanocoenia 23.2375 -81.046667 24.768957 0.000263
Undaria 23.2375 -81.046667 24.768957 0.001044
Xestospongia 23.2375 -81.046667 24.768957 0.000000

rerddap’s info request does not have enough metadata about the variables to explain the blank, and most abundant, genus. Checking the sever did not help figure that out. We’ll remove that for now to deal with only those that are identified.

cremp_1996 = cremp_1996.loc[cremp_1996["genus"] != ""]

There are also many genus with zero biomass count. In this example we’ll choose to do a biased analysis of occurrence and eliminate those where nothing was observed.

# Filter zero values (nothing was observed).

cremp_1996 = cremp_1996.loc[cremp_1996["organismQuantity"] > 0]

Now we can check the quantificationValue average by genus.

%matplotlib inline

avg = cremp_1996.groupby("genus").mean()  # re-compute the "biased" average.
ax = avg["organismQuantity"].plot(kind="bar")
../../../_images/2017-11-30-rerddap_17_0.png

and habitat.

ax = cremp_1996.groupby("habitat").mean()["organismQuantity"].plot(kind="bar")
../../../_images/2017-11-30-rerddap_19_0.png

It seems the most of the biomass was found around the Dendrogyra genus in Patch Reef habitats. But where are those Coral Reefs? How is the distribution of top three species with more biomass around them? With a pandas DataFrame it is easy to group the data by location and count the genus occurrence based on it.

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter


def make_plot():
    bbox = [-82, -80, 24, 26]
    projection = ccrs.PlateCarree()

    fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=dict(projection=projection))

    ax.set_extent(bbox)

    land = cfeature.NaturalEarthFeature(
        "physical", "land", "10m", edgecolor="face", facecolor=[0.85] * 3
    )

    ax.add_feature(land, zorder=0)
    ax.coastlines("10m", zorder=1)

    ax.set_xticks(np.linspace(bbox[0], bbox[1], 3), crs=projection)
    ax.set_yticks(np.linspace(bbox[2], bbox[3], 3), crs=projection)
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    return fig, ax
count = (
    cremp_1996.loc[cremp_1996["genus"] == "Acropora"]
    .groupby(["longitude", "latitude"])
    .count()
    .reset_index()
)


fig, ax = make_plot()
c = ax.scatter(
    count["longitude"],
    count["latitude"],
    s=200,
    c=count["genus"],
    alpha=0.5,
    cmap=plt.cm.get_cmap("viridis_r", 6),
    zorder=3,
)
cbar = fig.colorbar(c, shrink=0.75, extend="both")
cbar.ax.set_ylabel("Genus occurrence count")
ax.set_title("Acropora");
../../../_images/2017-11-30-rerddap_22_0.png
count = (
    cremp_1996.loc[cremp_1996["genus"] == "Dendrogyra"]
    .groupby(["longitude", "latitude"])
    .count()
    .reset_index()
)


fig, ax = make_plot()
c = ax.scatter(
    count["longitude"],
    count["latitude"],
    s=200,
    c=count["genus"],
    alpha=0.5,
    cmap=plt.cm.get_cmap("viridis_r", 6),
    zorder=3,
)
cbar = fig.colorbar(c, shrink=0.75, extend="both")
cbar.ax.set_ylabel("Genus occurrence count")
ax.set_title("Dendrogyra");
../../../_images/2017-11-30-rerddap_23_0.png
count = (
    cremp_1996.loc[cremp_1996["genus"] == "Orbicella"]
    .groupby(["longitude", "latitude"])
    .count()
    .reset_index()
)


fig, ax = make_plot()
c = ax.scatter(
    count["longitude"],
    count["latitude"],
    s=200,
    c=count["genus"],
    alpha=0.5,
    cmap=plt.cm.get_cmap("viridis_r", 6),
    zorder=3,
)
cbar = fig.colorbar(c, shrink=0.75, extend="both")
cbar.ax.set_ylabel("Genus occurrence count")
ax.set_title("Orbicella");
../../../_images/2017-11-30-rerddap_24_0.png

This demonstration showed the power of mixing Python and R to reduce developer time and allow the research to focus on the data and not the programming language.