{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d6f5c88",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright 2025 Google LLC\n",
    "#\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "#     https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccd3b043",
   "metadata": {},
   "source": [
    "# Data Engineering with GeoSpacial data\n",
    "\n",
    "![Cover](../../../docs/images/assessing_risks_geospatial_bqml/cover.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "988b9450",
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "source": [
    "| Author |\n",
    "| --- |\n",
    "| [Luis Gerardo Baeza](https://github.com/lgbaeza) |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "RwVIxeI9kIP2",
   "metadata": {
    "id": "RwVIxeI9kIP2"
   },
   "source": [
    "## Scenario\n",
    "For agribusinesses and their financial investors and insurers, accurately quantifying the financial risk from a production shortage is a big analytical challenge and an important one to solve. The direct impact of a fire or a drought is on the crop yield potential loss. These events lead to revenue loss, disrupt supply chain operations, and have financial consequences.\n",
    "\n",
    "## Goal\n",
    "On this notebook, we use the vast repository of public satellite imagery and environmental data available through Google Earth Engine, which is available as rasterized data directly within BigQuery. As an example, we will walk you through analyzing the wine production in Chile, to create dynamic wildfire risk intelligence system. By joining this environmental data with the geographic locations of vineyards and other assets, stakeholders can generate precise risk assessments for specific areas of interest.\n",
    "\n",
    "## High-level flow\n",
    "* **Access Chilean Regions of Interest:** Use Overture Maps in BigQuery Sharing to access regional data. Combine this data with official wine production statistics to map regions of interest and their geographical boundaries using the GEOGRAPHY data type in BigQuery.\n",
    "* **Explore Environmental Data:** Quickly explore environmental characteristics using Earth Engine’s public datasets. Analyze elevation, wildfires history, water irrigation, and availability, as well as Net Primary Production using datasets such as the NASA Digital Elevation Model, GFSAD Cropland, and GLOBathy Global Lakes Bathymetry.\n",
    "* **Perform Raster Analysis with ST_RegionStats():** Apply the ST_RegionStats() function to compute aggregate statistics (e.g., mean, sum, standard deviation) for raster data within Chilean regions defined by BigQuery geographies. This allows for summarizing raster insights over vector boundaries.\n",
    "* **Visualize and Interpret Results:** Export or connect BigQuery results to geospatial visualization tools (e.g., GeoViz, Data Studio, Looker, Carto) to interpret and communicate findings effectively.\n",
    "* **Get actionable insights:** Develop a Machine Learning model to predict risk areas and take timely actions.\n",
    "\n",
    "## Going further\n",
    "This is just a case study that might get your attention if you are a wine enthusiast, however, the same technology stack can solve complex location-based problems in a wide variety of companies and industries, from logistics network optimization to mineral exploration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ss_T0l3FkLf5",
   "metadata": {
    "id": "ss_T0l3FkLf5"
   },
   "source": [
    "## Step 0: Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ee0475f",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install bigframes pandas matplotlib google-cloud-bigquery -q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "510e2332",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext google.cloud.bigquery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2909e92",
   "metadata": {},
   "outputs": [],
   "source": [
    "!gcloud services enable earthengine.googleapis.com"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "DuJ9BovypV5s",
   "metadata": {
    "id": "DuJ9BovypV5s"
   },
   "outputs": [],
   "source": [
    "from google.cloud import bigquery\n",
    "import pandas as pd\n",
    "import bigframes.pandas as bpd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely import wkt\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "B2GLhIjeqtPr",
   "metadata": {
    "id": "B2GLhIjeqtPr"
   },
   "source": [
    "Configure a Google Cloud Project ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "z_qQ1L8pQj5F",
   "metadata": {
    "id": "z_qQ1L8pQj5F"
   },
   "outputs": [],
   "source": [
    "GOOGLE_CLOUD_PROJECT = \"\" # @param {type:\"string\"}\n",
    "if not GOOGLE_CLOUD_PROJECT: \n",
    "    GOOGLE_CLOUD_PROJECT = input(\"Enter project id\")\n",
    "\n",
    "BIGQUERY_DATASET = 'geospacial_data' # @param {type:\"string\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "KVj3bJjNSeBv",
   "metadata": {
    "id": "KVj3bJjNSeBv"
   },
   "outputs": [],
   "source": [
    "bq_client = bigquery.Client( project = GOOGLE_CLOUD_PROJECT )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9GTBaeCUusOc",
   "metadata": {
    "id": "9GTBaeCUusOc"
   },
   "outputs": [],
   "source": [
    "def quick_plot_map( df, geo_col, columns_color, size_x = 8, size_y = 8, subplots_x = 1, subplots_y = 1, title = None ):\n",
    "  gdf = gpd.GeoDataFrame( df, geometry = geo_col)\n",
    "  fig, ax = plt.subplots( subplots_y, subplots_x, figsize = ( size_y, size_x ) )\n",
    "\n",
    "  if subplots_x * subplots_y == 1:\n",
    "      ax_flat = [ ax ]\n",
    "  else:\n",
    "    ax_flat = ax.flatten( )\n",
    "\n",
    "  for i, col_color in enumerate( columns_color ):\n",
    "    current_ax = ax_flat[ i ]\n",
    "    gdf.plot(\n",
    "        column = col_color,\n",
    "        ax = current_ax,\n",
    "        legend = True,\n",
    "        cmap = 'inferno',\n",
    "        legend_kwds = { 'label' : col_color, 'orientation': \"horizontal\" }\n",
    "    )\n",
    "\n",
    "    # title and labels\n",
    "    if title is None:\n",
    "      current_ax.set_title( col_color.upper( ), fontdict = { 'fontsize': '12', 'fontweight': '3' } )\n",
    "    else:\n",
    "      current_ax.set_title( title, fontdict = { 'fontsize': '12', 'fontweight': '3' } )\n",
    "\n",
    "    current_ax.set_xlabel( 'Longitude' )\n",
    "    current_ax.set_ylabel( 'Latitude' )\n",
    "\n",
    "  # Hide unused subplots\n",
    "  for j in range(len(columns_color), len(ax_flat)):\n",
    "    ax_flat[j].set_visible(False)\n",
    "\n",
    "  plt.show( )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "uAdNRAmXSI7J",
   "metadata": {
    "id": "uAdNRAmXSI7J"
   },
   "source": [
    "## Step 1: Load data\n",
    "For the demonstration, we will use both publicly available data and a private-sourced table with the wine production in Chile by Region:\n",
    "- a) **Geographical regions in Chile** from the overture_maps dataset available in bigquery-public-data.\n",
    "- b) **Raster data** in Earth Engine (satellite images) from where we will extract  earth characteristics in regions of interest in Chile, like elevation, water irrigation and inland water body depth.\n",
    "- c) **Net Primary Production** (NPP) available in modis_terra_net_primary_production dataset of bigquery-public-data: a global measure of vegetation productivity.\n",
    "- d) **Wine production in Chile by region** within in a custom provided CSV file.\n",
    "- e) **Fire risk for Chilean regions** available in the nasa_wildfire.past_week table of bigquery-public-data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "kexKpwNtmvjr",
   "metadata": {
    "id": "kexKpwNtmvjr"
   },
   "source": [
    "Creation of a BigQuery dataset for the sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "Q0Kmy6fLmige",
   "metadata": {
    "id": "Q0Kmy6fLmige"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  CREATE SCHEMA IF NOT EXISTS { BIGQUERY_DATASET };\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "-fCyewJNmqdH",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "-fCyewJNmqdH",
    "outputId": "c8d2712e-4da4-4aa9-86d2-f8662de4e99a"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sIsAJ_HpnJDe",
   "metadata": {
    "id": "sIsAJ_HpnJDe"
   },
   "source": [
    "**a) Regions in Chile**\n",
    "\n",
    "The [overture maps dataset](https://docs.overturemaps.org/) in BigQuery public data provides free and open map data with representations of human settlements, such as countries, regions, states, cities, and even neighborhoods. The *geometry* column is a BigQuery Geography data type containing polygons with the boundaries for political divisions that we need, in this case, Chilean regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8YFetDG9nHpE",
   "metadata": {
    "id": "8YFetDG9nHpE"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "CREATE OR REPLACE TABLE `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.cl_regions` AS\n",
    "SELECT\n",
    "    names.primary AS region_name\n",
    "  , region\n",
    "  , st_Union_agg( geometry ) AS geometry\n",
    "FROM\n",
    "  `bigquery-public-data.overture_maps.division_area`\n",
    "WHERE\n",
    "  country = 'CL' AND subtype = 'region'\n",
    "  GROUP BY 1,2\n",
    "  ORDER BY region\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "w2Lq77BwnoU5",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "w2Lq77BwnoU5",
    "outputId": "bb9fd378-1513-442c-fd0b-f922a00ae5aa"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ttcokWasxHFw",
   "metadata": {
    "id": "ttcokWasxHFw"
   },
   "source": [
    "**b) Earth features from raster data in Earth Engine**\n",
    "\n",
    "Using Earth Engine's public datasets, we can acquire elevation, wildfires, water irrigation, and availability insights.\n",
    "* The [**Digital Elevation Model**](https://developers.google.com/earth-engine/datasets/catalog/USGS_SRTMGL1_003#description) provides valuable elevation data, which is relevant for wine grape growing..\n",
    "* The [**Cropland dataset**](https://developers.google.com/earth-engine/datasets/catalog/USGS_GFSAD1000_V1) provides the *landcover* band with information about water irrigation used for crops. A number between 1 - 2 means that either major or minor irrigation is performed. While a 0 value means there are no croplands. Water irrigation might influence water availability for wine acvitiy.\n",
    "* The [**Bathymetry**](https://developers.google.com/earth-engine/datasets/catalog/projects_sat-io_open-datasets_GLOBathy_GLOBathy_bathymetry) dataset contains detailed bathymetric maps by integrating maximum depth estimates comprising data on over 1.4 million waterbodies globally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mpisYG7qxTSS",
   "metadata": {
    "id": "mpisYG7qxTSS"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  CREATE OR REPLACE TABLE `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.earth_features_cl` AS\n",
    "  SELECT\n",
    "    region\n",
    "    , region_name\n",
    "    , ST_RegionStats(\n",
    "          geometry\n",
    "        , 'ee://USGS/SRTMGL1_003'                                   -- Digital elevation model (DEM)\n",
    "        , 'elevation'\n",
    "      ) AS elevation\n",
    "    , ST_RegionStats(\n",
    "          geometry\n",
    "        , 'USGS/GFSAD1000_V1'                                      -- Cropland\n",
    "        , 'landcover'\n",
    "    ) AS landcover\n",
    "    , ST_RegionStats(\n",
    "          geometry\n",
    "        , 'projects/sat-io/open-datasets/GLOBathy/GLOBathy_bathymetry'     -- Bathymetry\n",
    "        , 'b1'\n",
    "    ) AS water\n",
    "  FROM\n",
    "    `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.cl_regions`;\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "TTjhSLwTycj0",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "TTjhSLwTycj0",
    "outputId": "9627c212-209b-4c77-ee4b-92fa13a15c7a"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "-sF4bQ68PSzG",
   "metadata": {
    "id": "-sF4bQ68PSzG"
   },
   "source": [
    "**c) Net Primary Production**\n",
    "\n",
    "The [**modis_terra_net_primary_production**](https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD17A3HGF) provides the annual Gross Primary Productivity (GPP) and Net Primary Production (NPP), fundamental measures in ecology. It represents the  amount of carbon assimilated by plants and other photosynthetic organisms. The GPP is the total amount of carbon fixed by plants through photosynthesis.  NPP is what's left after GPP is used for the plant's own energy needs (respiration).\n",
    "\n",
    "NPP is relevant for wine growing because it indicates the overall health and vigor of the grapevine, as it measures the total amount of carbon the plant has converted into biomass.\n",
    "\n",
    "We will use the data already exported from Earth Engine into BigQuery, process that can be achieved with [Export.table.toBigQuery](https://developers.google.com/earth-engine/guides/exporting_to_bigquery)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "EMBR-iCNPVF-",
   "metadata": {
    "id": "EMBR-iCNPVF-"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  CREATE OR REPLACE TABLE `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.npp` AS\n",
    "    SELECT\n",
    "        region_name\n",
    "      , region\n",
    "      , year\n",
    "      , AVG( npp.npp ) AS avg_npp\n",
    "    FROM\n",
    "      `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.cl_regions`\n",
    "    JOIN\n",
    "      `bigquery-public-data.modis_terra_net_primary_production.MODIS_MOD17A3HGF` AS npp\n",
    "    ON ST_WITHIN( npp.geography, geometry )\n",
    "    GROUP BY 1, 2, 3\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "kgECbbV3PPjh",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "kgECbbV3PPjh",
    "outputId": "f9e9d80c-2fed-4368-d4f3-c5ef61d83765"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "jK0ZAYduioq1",
   "metadata": {
    "id": "jK0ZAYduioq1"
   },
   "source": [
    "**d) Wine production in Chile**\n",
    "\n",
    "The wine production is a dataset compiled from the official statistics published by the [Chilean authority SAG](https://www.sag.gob.cl) for years 2022, 2023, and 2024."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "onOmhGlwQVBT",
   "metadata": {
    "id": "onOmhGlwQVBT"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "-- Create the table\n",
    "CREATE OR REPLACE TABLE `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.cl_wine_production` (\n",
    "  region_name STRING\n",
    "  , region STRING\n",
    "  , liters_2024 INT64\n",
    "  , liters_2023 INT64\n",
    "  , liters_2022 INT64\n",
    ");\n",
    "\n",
    "-- Insert rows into the table\n",
    "INSERT INTO `{ GOOGLE_CLOUD_PROJECT }.{ BIGQUERY_DATASET }.cl_wine_production`\n",
    "  (region_name, region, liters_2024, liters_2023, liters_2022)\n",
    "VALUES\n",
    "  ('Tarapaca', 'CL-TA', 6779, 3434, 6923),\n",
    "  ('Antofagasta', 'CL-AN', 13897, 0, 0),\n",
    "  ('Atacama', 'CL-AT', 107494, 218137, 37820),\n",
    "  ('Coquimbo', 'CL-CO', 30816166, 72994218, 55611096),\n",
    "  ('Valparaiso', 'CL-VS', 19455875, 21032596, 19558904),\n",
    "  ('Metropolitana', 'CL-RM', 76852141, 87764852, 102801503),\n",
    "  ('Lib. Bernardo O Higgins', 'CL-LI', 327995699, 363757625, 420107504),\n",
    "  ('Maule', 'CL-ML', 464465227, 543907293, 624385653),\n",
    "  ('Nuble', 'CL-NB', 10647390, 13335017, 21513154),\n",
    "  ('Bio-Bio', 'CL-BI', 289000, 0, 332301),\n",
    "  ('Araucania', 'CL-AR', 9125, 10490, 9270),\n",
    "  ('Los Rios', 'CL-LR', 2025, 6669, 5853),\n",
    "  ('Los Lagos', 'CL-LL', 2005, 0, 0),\n",
    "  ('Aysen', 'CL-AI', 950, 1147, 0),\n",
    "  ('Magallanes y de la Antartica Chilena', 'CL-MA', 0, 0, 0),\n",
    "  ('Arica y Parinacota', 'CL-AP', 0, 0, 0)\n",
    ";\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eu4JH30N0RC0",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "eu4JH30N0RC0",
    "outputId": "12f6ea2f-6e88-4757-e880-4316c771978b"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "KCFhYzQTuz3C",
   "metadata": {
    "id": "KCFhYzQTuz3C"
   },
   "source": [
    "**e) Fire risk for Chilean regions**\n",
    "\n",
    "Using the *bigquery-public-data.nasa_wildfire.past_week* we compute a relative fire risk index for the Chilean regions.\n",
    "\n",
    "The dataset contains the bright_ti4 and bright_ti5 fields, which correspond to the brightness temperature measurements from the VIIRS (Visible Infrared Imaging Radiometer Suite) sensor onboard NASA satellites. The key difference lies in the wavelength each channel measures. The difference between bright_ti4 and bright_ti5 is a strong indicator of a real fire.\n",
    "\n",
    "We calculate a per-fire intensity score incorporating both brightness channels and also the FRP that quantifies the rate at which a fire releases energy in the form of electromagnetic radiation.\n",
    "\n",
    "Notice that we are using the BigQuery Geography function [ST_INTERSECTS](https://cloud.google.com/bigquery/docs/reference/standard-sql/geography_functions#st_intersects) to match each fire recorded in the *nasa_wildfire.past_week* table with the chilean regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tZDKj6FGtklp",
   "metadata": {
    "id": "tZDKj6FGtklp"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  CREATE OR REPLACE TABLE `{ BIGQUERY_DATASET }.fire_risk_cl` AS\n",
    "    WITH fire_analysis AS (\n",
    "      WITH fires AS (\n",
    "            SELECT\n",
    "              ST_GEOGPOINT(longitude, latitude) AS longlat\n",
    "              , *\n",
    "            FROM `bigquery-public-data.nasa_wildfire.past_week`\n",
    "          )\n",
    "          SELECT\n",
    "              region\n",
    "            , count(1) fire_count\n",
    "            , AVG ( (fires.bright_ti4 - fires.bright_ti5) * 0.2 + fires.frp * 0.8 )\n",
    "              fire_intensity\n",
    "          FROM { BIGQUERY_DATASET }.cl_regions regions\n",
    "          JOIN fires\n",
    "            ON ST_INTERSECTS(fires.longlat, geometry)\n",
    "            GROUP BY region\n",
    "    )\n",
    "        SELECT\n",
    "            regions.region\n",
    "          , regions.region_name\n",
    "          , IF(fire_count is NULL, 0, fire_count) fire_count\n",
    "          , IF(fire_intensity is NULL, 0, fire_intensity) fire_intensity\n",
    "          , IF(fire_count IS NULL, 0, 100 * PERCENT_RANK() OVER (ORDER BY fire_intensity * fire_count)) fire_index\n",
    "        FROM fire_analysis\n",
    "        RIGHT JOIN { BIGQUERY_DATASET }.cl_regions regions\n",
    "        ON fire_analysis.region = regions.region\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2BuiKLBPtvfW",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "2BuiKLBPtvfW",
    "outputId": "a6a0cebe-cdc4-42d9-bd5f-19a7df21ac5d"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "F0PElEtaGuv0",
   "metadata": {
    "id": "F0PElEtaGuv0"
   },
   "source": [
    "## Step 2: Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4lTzqQxoHT4V",
   "metadata": {
    "id": "4lTzqQxoHT4V"
   },
   "source": [
    "Now, a table is prepared with all the features we wanto to include in the analysis of the potential correlations and risks to the wine production, including wild fire risk, earth features. This features are put together along the Geography field to be able to analyze it and visualize them interactively in a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "_6QZwycQHBq9",
   "metadata": {
    "id": "_6QZwycQHBq9"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "CREATE OR REPLACE TABLE { BIGQUERY_DATASET }.wine_data_analysis AS\n",
    "SELECT\n",
    "    cl_regions.*\n",
    "  , fire_risk_cl.fire_index\n",
    "  , npp.avg_npp\n",
    "  , earth_features_cl.elevation.mean elevation_mean\n",
    "  , earth_features_cl.landcover.mean landcover_mean\n",
    "  , earth_features_cl.water.mean water_mean\n",
    "  , cl_wine_production.liters_2024\n",
    "FROM { BIGQUERY_DATASET }.cl_regions\n",
    "INNER JOIN { BIGQUERY_DATASET }.fire_risk_cl USING (region)\n",
    "INNER JOIN { BIGQUERY_DATASET }.npp USING (region)\n",
    "INNER JOIN { BIGQUERY_DATASET }.earth_features_cl USING (region)\n",
    "INNER JOIN { BIGQUERY_DATASET }.cl_wine_production USING (region)\n",
    "WHERE npp.year = (SELECT max(npp.year) FROM { BIGQUERY_DATASET }.npp)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "y800m7_JJEVx",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "y800m7_JJEVx",
    "outputId": "75d2975f-297c-4f54-aacc-f68e03852e7f"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "v2d_23JjoQQI",
   "metadata": {
    "id": "v2d_23JjoQQI"
   },
   "source": [
    "## Step 3: Exploring the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "N2bQXajdEe_5",
   "metadata": {
    "id": "N2bQXajdEe_5"
   },
   "source": [
    "Exploring the wine production using the Geography field with the boundaries of our interests (Regions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9PDIeH0AoLwI",
   "metadata": {
    "id": "9PDIeH0AoLwI"
   },
   "outputs": [],
   "source": [
    "SQL = f\"SELECT region_name, geometry, liters_2024 FROM { BIGQUERY_DATASET }.wine_data_analysis ORDER BY liters_2024 DESC LIMIT 10\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1j4HWoUdoPfZ",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 426
    },
    "id": "1j4HWoUdoPfZ",
    "outputId": "bba89746-9665-4630-cd58-b9276f0d064a"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "s-VXCVDdUOcP",
   "metadata": {
    "id": "s-VXCVDdUOcP"
   },
   "source": [
    " We use [BigQuery Geo Viz](https://bigquerygeoviz.appspot.com/) to visualize it on a map. We can effectively see that \"Maule\" region is the greatest wine productor in 2024 with 464,465,227 liters.  You can use the output of the following cell as the input query if you want to try it on your own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "AWlvx5vepmrl",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "id": "AWlvx5vepmrl",
    "outputId": "40e4912c-b5b6-448a-f844-56a4135aab9f"
   },
   "outputs": [],
   "source": [
    "f\"\"\"\n",
    "  SELECT region_name, geometry, liters_2024\n",
    "  FROM { BIGQUERY_DATASET }.wine_data_analysis\n",
    "\"\"\".replace(\"\\n\", \"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ivxeT5suLre",
   "metadata": {
    "id": "4ivxeT5suLre"
   },
   "source": [
    "\n",
    "\n",
    " To reproduce this map, you will need to setup the styles in the following way:\n",
    "* **fillColor** set to Data-driven with a linear fill using the liters_2024 field. The Domain must be set to 0 and the maximum, in this case, 464465227. The range must be set to #fff and #ab0000.\n",
    "* **fillOpacity** set to Data-driven with a linear fill using the liters_2024. The Domain must be set to 0 and the maximum, in this case, 464465227. The range must be set to 0 and 1.\n",
    "![BigQuery GeoViz](../../../docs/images/assessing_risks_geospatial_bqml/bq-regions-wine-product.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a-FPjIl2x1kW",
   "metadata": {
    "id": "a-FPjIl2x1kW"
   },
   "source": [
    "As an alternative, you can visualize it directly in Colab using the open source package [Geopandas](https://geopandas.org/en/stable/). For that, we can access the BigQuery data with using **BigFrames.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "_4eq-BGEpTme",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 287
    },
    "id": "_4eq-BGEpTme",
    "outputId": "74e83a88-e1f6-4627-ecd5-aa4280306436"
   },
   "outputs": [],
   "source": [
    "wine_production_by_region = bpd.read_gbq( f\"\"\"\n",
    "  SELECT region, geometry, liters_2024\n",
    "  FROM { BIGQUERY_DATASET }.wine_data_analysis\n",
    "\"\"\")\n",
    "wine_production_by_region.head( )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "Onn-pbwUSjdl",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 657
    },
    "id": "Onn-pbwUSjdl",
    "outputId": "afb12fb5-998b-425a-aa3e-b4a89546a0fb"
   },
   "outputs": [],
   "source": [
    "quick_plot_map(\n",
    "    wine_production_by_region.to_pandas()\n",
    "    , columns_color = [ 'liters_2024' ]\n",
    "    , geo_col = 'geometry'\n",
    "    , title = 'Wine production per Region in 2024'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hFUtkltdxC1D",
   "metadata": {
    "id": "hFUtkltdxC1D"
   },
   "source": [
    "Let's also visualize the earth features that we acquire from the raster data and the bigquery-public-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "Ot-Ogu6TWFAG",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "id": "Ot-Ogu6TWFAG",
    "outputId": "7f63d9b2-0084-4308-83e4-6385c573ec2b"
   },
   "outputs": [],
   "source": [
    "chile_features = bpd.read_gbq( f\"\"\"\n",
    "  SELECT\n",
    "    region_name\n",
    "    , geometry\n",
    "    , elevation_mean\n",
    "    , water_mean\n",
    "    , avg_npp\n",
    "    , landcover_mean\n",
    "    , fire_index\n",
    "  FROM { BIGQUERY_DATASET }.wine_data_analysis\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mC2Xrnrbmmb6",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 621
    },
    "id": "mC2Xrnrbmmb6",
    "outputId": "9aca2327-c2a4-40a2-b4da-2ce8ca0cc94c"
   },
   "outputs": [],
   "source": [
    "columns_to_plot = [\"elevation_mean\", \"water_mean\", \"avg_npp\", \"landcover_mean\", \"fire_index\"]\n",
    "quick_plot_map(\n",
    "    chile_features.to_pandas()\n",
    "    , columns_color = columns_to_plot\n",
    "    , geo_col = 'geometry'\n",
    "    , subplots_x = 3\n",
    "    , subplots_y = 2\n",
    "    , size_x = 7\n",
    "    , size_y = 8\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mCE7R-hHfiUq",
   "metadata": {
    "id": "mCE7R-hHfiUq"
   },
   "source": [
    "## Step 4: Unlocking Hidden insights with a K-Means model using BigQuery ML"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ETXzC84ngLAD",
   "metadata": {
    "id": "ETXzC84ngLAD"
   },
   "source": [
    "Now we will train a K-Means clustering model as an example method to segment land holdings into distinct groups based on the environmental characteristics already acquired from Earth Engine and BigQuery public data. As an unsupervised learning model, it automatically identifies natural groupings of regions without needing a pre-defined label. A company can use the resulting clusters to discover, for example, a \"drought-prone\" cluster that requires specific irrigation strategies, a \"high-yield\" cluster that warrants investment, or a \"fire-vulnerable\" cluster that needs targeted prevention programs. By moving beyond simple regional boundaries, this model provides a nuanced and data-driven perspective on land management, allowing for the creation of tailored, efficient, and highly effective agricultural strategies.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "TMvfx_wu6gFc",
   "metadata": {
    "id": "TMvfx_wu6gFc"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  CREATE OR REPLACE MODEL\n",
    "    `{ BIGQUERY_DATASET }.kmeans_model`\n",
    "  OPTIONS(\n",
    "    model_type = 'KMEANS',\n",
    "    num_clusters = 4,\n",
    "    distance_type = 'EUCLIDEAN'\n",
    "  ) AS\n",
    "  SELECT\n",
    "    elevation_mean\n",
    "    , avg_npp\n",
    "    , landcover_mean\n",
    "    , water_mean\n",
    "    , fire_index\n",
    "    , liters_2024\n",
    "    , ST_CENTROID(geometry) AS geometry_point\n",
    "  FROM\n",
    "    `{ BIGQUERY_DATASET }.wine_data_analysis`;\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98mjETJZfhIC",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "98mjETJZfhIC",
    "outputId": "cc46c4d4-3050-46d4-df75-6705bfc58f98"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fpAQAG30kzl5",
   "metadata": {
    "id": "fpAQAG30kzl5"
   },
   "source": [
    "After aproximately one minute, we will find that the model was trained and we can see under **Execution details** some information about the iterations and duration.\n",
    "![bqml_model_trained](../../../docs/images/assessing_risks_geospatial_bqml/bqml_model_trained.png)\n",
    "\n",
    "Under Resutls, if you choose **Go To Model**, you can also find stats about the clusters identified.\n",
    "![bqml_cluster_stats](../../../docs/images/assessing_risks_geospatial_bqml/bqml_cluster_stats.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "yxapNxR8ohgF",
   "metadata": {
    "id": "yxapNxR8ohgF"
   },
   "source": [
    "## Step 5: Actionable cluster insights for strategic planning\n",
    "This analysis identifies four distinct operational profiles for the wine producing regions, offering clear guidance for resource allocation and strategic planning.  The highest-producing clusters also carry the greatest risk.\n",
    "\n",
    "* **Cluster 3** is your most profitable group, with the highest annual production of over 18 billion liters, but it also has a significant fire index and high elevation, which indicates a need for specialized terrain management and robust risk mitigation.\n",
    "* **Cluster 1** is also a high-yield, high-risk group, producing over 267 million liters, making it a priority for fire prevention and resource protection.\n",
    "* **Cluster 2** presents a low-risk, high-potential opportunity for future business development, due to its ample water and zero fire activity.\n",
    "* Conversely, **Cluster 4** is an anomaly and an area for immediate investigation: despite having the highest average water availability, it has very low production. This indicates a potential inefficiency or different land use that could be optimized to increase output, representing a significant opportunity for growth.\n",
    "\n",
    "Let's find below which regions belong to the clusters using the *ML.PREDICT* function to trigger the model inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "s5NCU_epfiqC",
   "metadata": {
    "id": "s5NCU_epfiqC"
   },
   "outputs": [],
   "source": [
    "SQL = f\"\"\"\n",
    "  WITH cluster_analysis as (\n",
    "    SELECT\n",
    "      region, region_name, centroid_id CLUSTER_ID, rank() OVER(PARTITION BY centroid_id ORDER BY liters_2024 DESC) rank\n",
    "    FROM\n",
    "      ML.PREDICT(MODEL `{ BIGQUERY_DATASET }.kmeans_model`,\n",
    "        (\n",
    "          SELECT\n",
    "            region\n",
    "            , region_name\n",
    "            , elevation_mean\n",
    "            , avg_npp\n",
    "            , landcover_mean\n",
    "            , water_mean\n",
    "            , fire_index\n",
    "            , liters_2024\n",
    "            , ST_CENTROID(geometry) AS geometry_point\n",
    "          FROM\n",
    "            `{ BIGQUERY_DATASET }.wine_data_analysis`\n",
    "        )\n",
    "      )\n",
    "      ORDER BY\n",
    "        CASE centroid_id\n",
    "          WHEN 3 THEN 1\n",
    "          WHEN 1 THEN 2\n",
    "          WHEN 2 THEN 3\n",
    "          ELSE 4\n",
    "      END\n",
    "    )\n",
    "    SELECT region_name, CLUSTER_ID\n",
    "    FROM cluster_analysis\n",
    "    WHERE rank <= 3 AND CLUSTER_ID IN (1, 3);\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "g1FZcRv_rQ5H",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 301
    },
    "id": "g1FZcRv_rQ5H",
    "outputId": "3a1e753a-05c0-4247-f64a-3ce07c5421d3"
   },
   "outputs": [],
   "source": [
    "%%bigquery\n",
    "$SQL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "F-ZkIWDMtp1Z",
   "metadata": {
    "id": "F-ZkIWDMtp1Z"
   },
   "source": [
    "**Conclusion**\n",
    "\n",
    "This analysis demonstrates how to pinpoint areas with both high agricultural output and environmental risks, such as wildfires. This insight is crucial for implementing proactive, data-driven strategies that protect investments and mitigate potential losses.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "xHFhaET7nggH",
   "metadata": {
    "id": "xHFhaET7nggH"
   },
   "source": [
    "## Cleanup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "lMw2eQAdnhbR",
   "metadata": {
    "id": "lMw2eQAdnhbR"
   },
   "outputs": [],
   "source": [
    "dataset_ref = bq_client.dataset(BIGQUERY_DATASET)\n",
    "try:\n",
    "    # Use delete_contents=True to delete all tables in the dataset.\n",
    "    # Use not_found_ok=True to avoid an error if the dataset doesn't exist.\n",
    "    bq_client.delete_dataset( dataset_ref, delete_contents = True, not_found_ok = False)\n",
    "    print(f\"Dataset '{BIGQUERY_DATASET}' in project '{GOOGLE_CLOUD_PROJECT}' deleted successfully.\")\n",
    "except Exception as e:\n",
    "    print(f\"Error deleting dataset '{BIGQUERY_DATASET}': {e}\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "datascience_book",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "myenv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}