Covariate_sampling_102pts.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c6886c-57ea-4fec-9c1e-a6802d747f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################################################################################################################\n",
    "######################## This script extracts 55 covariates values at each of the soil algae data point position #########################\n",
    "##########################################################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "167bd8db-f5d9-4d4e-9ec0-60d1ef3f1bdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://accounts.google.com/o/oauth2/auth?client_id=517222506229-vsmmajv00ul0bs7p89v5m89qs8eb9359.apps.googleusercontent.com&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fearthengine+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fdevstorage.full_control&redirect_uri=urn%3Aietf%3Awg%3Aoauth%3A2.0%3Aoob&response_type=code&code_challenge=JiUs9rWOluOAe3Ykrzo0ZLf1U8ppM1QMvK3lPqapOKg&code_challenge_method=S256>https://accounts.google.com/o/oauth2/auth?client_id=517222506229-vsmmajv00ul0bs7p89v5m89qs8eb9359.apps.googleusercontent.com&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fearthengine+https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fdevstorage.full_control&redirect_uri=urn%3Aietf%3Awg%3Aoauth%3A2.0%3Aoob&response_type=code&code_challenge=JiUs9rWOluOAe3Ykrzo0ZLf1U8ppM1QMvK3lPqapOKg&code_challenge_method=S256</a></p>\n",
       "        <p>The authorization workflow will generate a code, which you\n",
       "        should paste in the box below</p>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Enter verification code:  4/1AX4XfWiTsXEThSEhjvC_7_XoAtq2gINsN5a2ATUKAPaHLVt_r4QllwxPXoQ\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Successfully saved authorization token.\n"
     ]
    }
   ],
   "source": [
    "import ee # Import the Python API package\n",
    "ee.Authenticate() # Authenticate to the Earth Engine servers\n",
    "ee.Initialize() # Initialize the API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a66e3f91-fba6-4ecd-b879-096ad6b69207",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geemap # The geemap Python package is built upon the ipyleaflet and folium packages and implements several methods for interacting with Earth Engine data layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "9acc2b45-6a01-467e-b12e-f989c3c3e282",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "10b6bada7bd74549839fba84be584357",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[0, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Toggl…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a map\n",
    "Map1 = geemap.Map(center=[0,0], zoom=2)\n",
    "Map1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "da7358d0-bcd3-478e-abe4-5efc0bef6951",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of the productivity data collection 102\n"
     ]
    }
   ],
   "source": [
    "# Loading sampling points stored in personnal repository\n",
    "ProductivityPts = ee.FeatureCollection('users/romainwalcker/Phototrophs_Productivity_AllPts')\n",
    "print('Size of the productivity data collection',ProductivityPts.size().getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7dcee8cc-cb98-497d-b1d8-136d20dde362",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add points to the map\n",
    "Map1.addLayer(ProductivityPts, {},'Sampling points')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "80e57468-cec8-48f3-af73-873297c59cf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading environmental variables (Table S3 of the paper)\n",
    "# Climate variables\n",
    "WorldClim = ee.Image('WORLDCLIM/V1/BIO');\n",
    "bio01 = WorldClim.select('bio01')# Annual Mean Temperature\n",
    "bio02 = WorldClim.select('bio02')# Mean Diurnal Range\n",
    "bio03 = WorldClim.select('bio03')# Isothermality\n",
    "bio04 = WorldClim.select('bio04')# Temperature Seasonality\n",
    "bio05 = WorldClim.select('bio05')# Max Temperature of Warmest Month\n",
    "bio06 = WorldClim.select('bio06')# Min Temperature of Coldest Month\n",
    "bio07 = WorldClim.select('bio07')# Temperature Annual Range (BIO5-BIO6)\n",
    "bio08 = WorldClim.select('bio08')# Mean Temperature of Wettest Quarter\n",
    "bio09 = WorldClim.select('bio09')# Mean Temperature of Driest Quarter\n",
    "bio10 = WorldClim.select('bio10')# Mean Temperature of Warmest Quarter\n",
    "bio11 = WorldClim.select('bio11')# Mean Temperature of Coldest Quarter\n",
    "bio12 = WorldClim.select('bio12')# Annual Precipitation\n",
    "bio13 = WorldClim.select('bio13')# Precipitation of Wettest Month\n",
    "bio14 = WorldClim.select('bio14')# Precipitation of Driest Month\n",
    "bio15 = WorldClim.select('bio15')# Precipitation Seasonality\n",
    "bio16 = WorldClim.select('bio16')# Precipitation of Wettest Quarter\n",
    "bio17 = WorldClim.select('bio17')# Precipitation of Driest Quarter\n",
    "bio18 = WorldClim.select('bio18')# Precipitation of Warmest Quarter\n",
    "bio19 = WorldClim.select('bio19')# Precipitation of Coldest Quarter\n",
    "\n",
    "TerraClim = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')\n",
    "aet = TerraClim.select('aet').mean()# Actual evapotranspiration\n",
    "tdef = TerraClim.select('def').mean()# Climate water deficit\n",
    "pdsi = TerraClim.select('pdsi').mean()# Palmer Drought Severity Index\n",
    "pet = TerraClim.select('pet').mean()# Reference evapotranspiration\n",
    "pr = TerraClim.select('pr').mean()# Precipitation accumulation\n",
    "ro = TerraClim.select('ro').mean()# Runoff \n",
    "srad = TerraClim.select('srad').mean()# Downward surface shortwave radiation\n",
    "swe = TerraClim.select('swe').mean()# Snow water equivalent\n",
    "tmmn = TerraClim.select('tmmn').mean()# Minimum temperature\n",
    "tmmx = TerraClim.select('tmmx').mean()# Maximum temperature\n",
    "vap = TerraClim.select('vap').mean()# Vapor pressure\n",
    "vpd = TerraClim.select('vpd').mean()# Vapor pressure deficit\n",
    "wspd = TerraClim.select('vs').mean().rename('wspd')# Wind-speed at 10m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a39862f4-99ea-4a36-be55-f31124b1adc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Soil variables\n",
    "USDA4A1H_M = ee.Image(\"OpenLandMap/SOL/SOL_BULKDENS-FINEEARTH_USDA-4A1H_M/v02\");\n",
    "bulk = USDA4A1H_M.select('b0').unmask().rename('bulk')# Soil bulk density at 0 cm depth\n",
    "USDA3A1A1A_M = ee.Image(\"OpenLandMap/SOL/SOL_CLAY-WFRACTION_USDA-3A1A1A_M/v02\");\n",
    "clay = USDA3A1A1A_M.select('b0').rename('clay')# Clay content at 0 cm depth\n",
    "USDA6A1C_M = ee.Image(\"OpenLandMap/SOL/SOL_ORGANIC-CARBON_USDA-6A1C_M/v02\");\n",
    "soc = USDA6A1C_M.select('b0').rename('soc')# Soil organic carbon content at 0 cm depth\n",
    "USDA4C1A2A_M = ee.Image(\"OpenLandMap/SOL/SOL_PH-H2O_USDA-4C1A2A_M/v02\");\n",
    "ph = USDA4C1A2A_M.select('b0').rename('ph')# Soil pH in H2O at 0 cm depth\n",
    "USDATT_M = ee.Image(\"OpenLandMap/SOL/SOL_TEXTURE-CLASS_USDA-TT_M/v02\");\n",
    "text = USDATT_M.select('b0').rename('text')# Soil texture classes\n",
    "USDA3A1A1A_M = ee.Image(\"OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02\");\n",
    "sand = USDA3A1A1A_M.select('b0').rename('sand')# Sand content\n",
    "USDA4B1C_M = ee.Image(\"OpenLandMap/SOL/SOL_WATERCONTENT-33KPA_USDA-4B1C_M/v01\");\n",
    "water = USDA4B1C_M.select('b0').rename('water')# Soil water content\n",
    "\n",
    "TerraClim = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE');\n",
    "moist = TerraClim.select('soil').mean().rename('moist')# Soil moisture, derived using a one-dimensional soil water balance model\n",
    "\n",
    "SMAP = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture');\n",
    "ssm = SMAP.select('ssm').mean()# Surface soil moisture\n",
    "susm = SMAP.select('susm').mean()# Subsurface soil moisture\n",
    "smp = SMAP.select('smp').mean()# Soil moisture profile\n",
    "ssma = SMAP.select('ssma').mean()# Surface soil moisture anomaly\n",
    "susma = SMAP.select('susma').mean()# Subsurface soil moisture anomaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8613850c-12fa-49a7-a073-a18f3944ca86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Topographic variables\n",
    "GMTED = ee.Image('USGS/GMTED2010')\n",
    "height = GMTED.select('be75').rename('height')# Elevation\n",
    "  \n",
    "Coordonates = ee.Image.pixelLonLat()\n",
    "lon = Coordonates.select('longitude').rename('lon')\n",
    "lat = Coordonates.select('latitude').rename('lat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1c6ee55e-d377-4abf-9b94-4cb737048f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vegetation variables\n",
    "MOD13Q1 = ee.ImageCollection('MODIS/006/MOD13Q1')\n",
    "ndvi = MOD13Q1.select('NDVI').mean().rename('ndvi')# Normalized Difference Vegetation Index\n",
    "evi = MOD13Q1.select('EVI').mean().rename('evi')# Enhanced Vegetation Index\n",
    "  \n",
    "MOD17A2H = ee.ImageCollection('MODIS/006/MOD17A2H')\n",
    "gpp = MOD17A2H.select('Gpp').mean().rename('gpp')# Gross primary production\n",
    "npp = MOD17A2H.select('PsnNet').mean().rename('npp')# Net photosynthesis\n",
    "  \n",
    "MCD15A3H = ee.ImageCollection('MODIS/006/MCD15A3H');\n",
    "lai = MCD15A3H.select('Lai').mean().rename('lai')# leaf area index\n",
    "\n",
    "biomes = ee.Image('OpenLandMap/PNV/PNV_BIOME-TYPE_BIOME00K_C/v01').rename('biomes')# Potential distribution of biomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "09ed6c78-78c1-45e4-8468-da6021358c29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Human variables\n",
    "GPW = ee.ImageCollection(\"CIESIN/GPWv411/GPW_Population_Density\").first()\n",
    "pop = GPW.select('population_density').rename('pop')# The estimated number of persons per square kilometer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "63e6e34a-37d0-41bf-bbfa-e227d3ba93a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stacking environmental variable in a multiband image\n",
    "compositeOfInterest = biomes.addBands(bio01).addBands(bio02).addBands(bio03).addBands(bio04).addBands(bio05).addBands(bio06).addBands(bio07).addBands(bio08).addBands(bio09).addBands(bio10).addBands(bio11).addBands(bio12).addBands(bio13).addBands(bio14).addBands(bio15).addBands(bio16).addBands(bio17).addBands(bio18).addBands(bio19).addBands(aet).addBands(tdef).addBands(pdsi).addBands(pet).addBands(pr).addBands(ro).addBands(srad).addBands(swe).addBands(tmmn).addBands(tmmx).addBands(vap).addBands(vpd).addBands(wspd).addBands(ndvi).addBands(evi).addBands(gpp).addBands(npp).addBands(lai).addBands(bulk).addBands(clay).addBands(soc).addBands(ph).addBands(text).addBands(sand).addBands(water).addBands(moist).addBands(ssm).addBands(susm).addBands(smp).addBands(ssma).addBands(susma).addBands(height).addBands(lat).addBands(lon).addBands(pop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "13943004-e44f-4fe5-aed0-d3be601eab38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of the sampled data collection 102\n",
      "First sample {'type': 'FeatureCollection', 'columns': {'Group': 'String', 'ID': 'Float', 'aet': 'Float<-32768.0, 32767.0>', 'bio01': 'Float', 'bio02': 'Float', 'bio03': 'Float', 'bio04': 'Float', 'bio05': 'Float', 'bio06': 'Float', 'bio07': 'Float', 'bio08': 'Float', 'bio09': 'Float', 'bio10': 'Float', 'bio11': 'Float', 'bio12': 'Float', 'bio13': 'Float', 'bio14': 'Float', 'bio15': 'Float', 'bio16': 'Float', 'bio17': 'Float', 'bio18': 'Float', 'bio19': 'Float', 'biomes': 'Short<0, 255>', 'bulk': 'Short<0, 255>', 'clay': 'Short<0, 255>', 'def': 'Float<-32768.0, 32767.0>', 'evi': 'Float<-32768.0, 32767.0>', 'gpp': 'Float<-32768.0, 32767.0>', 'height': 'Short', 'lai': 'Float<0.0, 255.0>', 'lat': 'Float', 'log_primpr': 'Float', 'lon': 'Float', 'moist': 'Float<-32768.0, 32767.0>', 'ndvi': 'Float<-32768.0, 32767.0>', 'npp': 'Float<-32768.0, 32767.0>', 'pdsi': 'Float<-32768.0, 32767.0>', 'pet': 'Float<-32768.0, 32767.0>', 'ph': 'Short<0, 255>', 'pop': 'Float', 'pr': 'Float<-32768.0, 32767.0>', 'primprod': 'Float', 'ro': 'Float<-32768.0, 32767.0>', 'sand': 'Short<0, 255>', 'smp': 'Float', 'soc': 'Short<0, 255>', 'srad': 'Float<-32768.0, 32767.0>', 'ssm': 'Float', 'ssma': 'Float', 'susm': 'Float', 'susma': 'Float', 'swe': 'Float<-32768.0, 32767.0>', 'system:index': 'String', 'text': 'Short<0, 255>', 'tmmn': 'Float<-32768.0, 32767.0>', 'tmmx': 'Float<-32768.0, 32767.0>', 'vap': 'Float<-32768.0, 32767.0>', 'vpd': 'Float<-32768.0, 32767.0>', 'water': 'Short<0, 255>', 'wspd': 'Float<-32768.0, 32767.0>'}, 'features': [{'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [11.99999832245454, 47.03000057856102]}, 'id': '00000000000000000002', 'properties': {'Group': 'Training', 'ID': 3, 'aet': 425.82142857142856, 'bio01': 16, 'bio02': 81, 'bio03': 31, 'bio04': 6075, 'bio05': 149, 'bio06': -106, 'bio07': 255, 'bio08': 92, 'bio09': -61, 'bio10': 92, 'bio11': -61, 'bio12': 1079, 'bio13': 148, 'bio14': 54, 'bio15': 35, 'bio16': 412, 'bio17': 171, 'bio18': 412, 'bio19': 171, 'biomes': 9, 'bulk': 100, 'clay': 17, 'def': 10.314814814814815, 'evi': 1849.7549019607843, 'gpp': 112.78500986193293, 'height': 2068, 'lai': 5.772375690607735, 'lat': 47.029381454500005, 'log_primpr': 3.735285827, 'lon': 12.001010946500003, 'moist': 154.07804232804233, 'ndvi': 3331.401960784314, 'npp': 92.33431952662721, 'pdsi': 3.1706349206349205, 'pet': 446.0224867724868, 'ph': 49, 'pop': 0, 'pr': 91.59126984126983, 'primprod': 42.900000003, 'ro': 48.992063492063494, 'sand': 41, 'smp': 0.8853220343589783, 'soc': 26, 'srad': 1392.228835978836, 'ssm': 22.223556518554688, 'ssma': 0.018492043018341064, 'susm': 37.97834014892578, 'susma': 0.024707498028874397, 'swe': 70.4973544973545, 'text': 7, 'tmmn': -9.404761904761905, 'tmmx': 81.87169312169313, 'vap': 637.9629629629629, 'vpd': 28.156084656084655, 'water': 38, 'wspd': 257.37566137566137}}]}\n"
     ]
    }
   ],
   "source": [
    "# Sampling\n",
    "# The code below allows for sampling the collection without dropping missing values for variables reduction\n",
    "ProductivityPtsSampled = compositeOfInterest.reduceRegions(reducer= ee.Reducer.first(),collection= ProductivityPts,scale= compositeOfInterest.select('soc').projection().nominalScale().getInfo())\n",
    "print('Size of the sampled data collection',ProductivityPtsSampled.size().getInfo())\n",
    "print('First sample', ProductivityPtsSampled.limit(1).getInfo());"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}