notebooks/R/full_export.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OBIS full export\n",
    "\n",
    "Sometimes you need to do analyses using big chunks of OBIS data. In those cases, downloading it through the API is not practical. OBIS also offer the full export of all OBIS occurrence data in `parquet` or `csv` formats.\n",
    "\n",
    "You can download the full export here (~15GB): https://obis.org/data/access/\n",
    "\n",
    "We will work with the `parquet` file. [Parquet](https://parquet.apache.org/) is a column-based performant and light format developed by Apache. We can work with `parquet` files on R and Python using the `arrow` [package.](https://arrow.apache.org/docs/r/index.html). You can learn more about parquet in [this tutorial](https://resources.obis.org/tutorials/arrow-obis/).\n",
    "\n",
    "Loading all this data at once in your computer is usually impossible, as it will not fit in memory. But that is not a problem! Arrow offer us many ways of loading just the structure of the file and perform filters without loading the dataset in memory. Let's explore how to do that.\n",
    "\n",
    "<span style=\"color: orange;\">NOTE: instead of downloading all the data, we will use the S3 file in this practical. However, it is an older copy. For the most up-to-date always download from the website.</span>\n",
    "\n",
    "You may need to install `arrow` and `h3jsr`. If that is the case, before calling `install.packages()` run those two lines:\n",
    "\n",
    "```\n",
    "install.packages(\"bspm\")\n",
    "suppressMessages(bspm::enable())\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(arrow)\n",
    "library(dplyr)\n",
    "library(ggplot2)\n",
    "\n",
    "ds_path <- \"s3://obis-products/exports/obis_20230726.parquet\"\n",
    "\n",
    "ds <- open_dataset(ds_path)\n",
    "\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call `arrow::open_dataset` to open a file (or a folder with multiple files) without loading it into memory. Note that you can just see the columns and data types.\n",
    "\n",
    "As a first example, we will select only the columns `scientificName`, `family`, and `date_year`. Then, we will filter by the family **Pomacentridae**. Finally, we will group by `scientificName` and `date_year` and count the number of records:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "pomac_records <- ds |> \n",
    "    select(scientificName, family, date_year) |> \n",
    "    filter(family == \"Pomacentridae\") |> \n",
    "    group_by(scientificName, date_year) |> \n",
    "    count() |> \n",
    "    collect()\n",
    "\n",
    "head(pomac_records)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now let's get the records for the \"Nemo\"s species, and plot in a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "sp_nemos <- robis::checklist(\"Amphiprion\")\n",
    "\n",
    "sp_nemos <- sp_nemos[sp_nemos$taxonRank == \"Species\",]\n",
    "\n",
    "sp_nemo_top4 <- sp_nemos[order(sp_nemos$records, decreasing = T), ][1:4,]\n",
    "\n",
    "nemo_data <- ds |> \n",
    "    select(scientificName, decimalLongitude, decimalLatitude) |> \n",
    "    filter(scientificName %in% sp_nemo_top4$scientificName) |> \n",
    "    collect()\n",
    "\n",
    "wrld <- rnaturalearth::ne_countries(returnclass = \"sf\")\n",
    "\n",
    "ggplot(nemo_data) +\n",
    "    geom_sf(data = wrld) +\n",
    "    geom_point(aes(x = decimalLongitude, y = decimalLatitude, color = scientificName)) +\n",
    "    facet_wrap(~ scientificName) +\n",
    "    theme_light()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although this works fine, for big queries it can be a bit slow, even if the file is local. Another approach that is very good is [`duckDB`](https://duckdb.org/docs/api/r.html). It not only enable to do the same queries, using SQL, but it also have a [spatial extension.](https://duckdb.org/docs/extensions/spatial/functions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(DBI)\n",
    "library(duckdb)\n",
    "library(glue)\n",
    "\n",
    "con <- dbConnect(duckdb())\n",
    "dbSendQuery(con, \"install spatial; load spatial;\")\n",
    "dbSendQuery(con, \"install httpfs; load httpfs;\")\n",
    "\n",
    "query <- glue(\n",
    "    \"\n",
    "    select scientificName, count(*) as number_records\n",
    "    from read_parquet('{ds_path}')\n",
    "    where family = 'Pomacentridae'\n",
    "    group by scientificName\n",
    "    \"\n",
    ")\n",
    "\n",
    "pom_fam <- dbGetQuery(con, query)\n",
    "\n",
    "pom_fam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for spatial queries is also straightforward. We just need to read the longitude and latitude as Spatial points and then use ST_Intersects to intersect with our geometry (the WKT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# And for the case of an area is much faster because it supports spatial queries!\n",
    "wkt_area <- \"POLYGON ((-79.189453 27.293689, -79.584961 23.765237, -75.9375 22.43134, -73.959961 24.726875, -74.750977 27.176469, -79.189453 27.293689))\"\n",
    "\n",
    "query <- glue::glue(\n",
    "    \"\n",
    "  select distinct scientificName\n",
    "  from read_parquet('{ds_path}')\n",
    "  where ST_Intersects(ST_Point(decimalLongitude, decimalLatitude), ST_GeomFromText('{wkt_area}'))\n",
    "  \"\n",
    ")\n",
    "\n",
    "all_sps_area <- dbGetQuery(con, query)\n",
    "\n",
    "head(all_sps_area)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}