Skip to content

Instantly share code, notes, and snippets.

@giswqs

giswqs/zonal statistics_tb.ipynb

Created September 22, 2020 21:38
Show Gist options
  • Save giswqs/1d7125d346642f21bb5eea9e4caa7838 to your computer and use it in GitHub Desktop.
Save giswqs/1d7125d346642f21bb5eea9e4caa7838 to your computer and use it in GitHub Desktop.
Download ZIP
zonal statistics
Display the source blob
Display the rendered blob
Raw
zonal statistics_tb.ipynb
{
"cells": [
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import ee\nimport geemap",
"execution_count": 26,
"outputs": []
},
{
"metadata": {
"scrolled": false,
"trusted": true
},
"cell_type": "code",
"source": "Map = geemap.Map()\nMap",
"execution_count": 27,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "Map(center=[40, -100], controls=(WidgetControl(options=['position'], widget=HBox(children=(ToggleButton(value=…",
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "d93cc0920389471093367a77f5b1e26b"
}
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Replace roi with your featureCollection"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "roi = ee.FeatureCollection('TIGER/2018/States') \\\n .filterBounds(ee.Geometry.Rectangle([-75.9311, 38.7880, -75.5467, 39.0956]))\nMap.centerObject(roi, 8)\nMap.addLayer(roi, {}, \"roi\")",
"execution_count": 28,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Composite image bands"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "NEDelevation = ee.Image('USGS/NED').clip(roi)\nMap.addLayer(NEDelevation, {'min':0, 'max':4000},'NED')",
"execution_count": 29,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "slope = ee.Terrain.slope(NEDelevation).clip(roi)\nMap.addLayer(slope,{'min':0,'max':45},'Slope')",
"execution_count": 30,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "aspect = ee.Terrain.aspect(NEDelevation).clip(roi)\nMap.addLayer(aspect,{'min':0,'max':360},'Aspect')",
"execution_count": 31,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def maskS2clouds(image):\n qa = image.select('QA60')\n cloudBitMask = 1 << 10\n cirrusBitMask = 1 << 11\n mask = qa.bitwiseAnd(cloudBitMask).eq(0) \\\n .And(qa.bitwiseAnd(cirrusBitMask).eq(0))\n return image.updateMask(mask).divide(10000)",
"execution_count": 32,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Create monthly (or daily, yearly) date sequence "
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# start_date = '2018-01-01'\n# end_date = '2018-06-30'\n# dates = geemap.date_sequence(start_date, end_date, 'month')",
"execution_count": 33,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "start_year = 2016\nend_year = 2020\nstart_month = 6\nend_month = 10\n\ndate_list = []\nfor year in range(start_year, end_year):\n for month in range(start_month, end_month):\n item = str(year) + '-' + str(month).zfill(2) + '-01'\n date_list.append(item)\nprint(date_list)",
"execution_count": 34,
"outputs": [
{
"output_type": "stream",
"text": "['2016-06-01', '2016-07-01', '2016-08-01', '2016-09-01', '2017-06-01', '2017-07-01', '2017-08-01', '2017-09-01', '2018-06-01', '2018-07-01', '2018-08-01', '2018-09-01', '2019-06-01', '2019-07-01', '2019-08-01', '2019-09-01']\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "dates = ee.List(date_list)",
"execution_count": 35,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "dates.getInfo()",
"execution_count": 36,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 36,
"data": {
"text/plain": "['2016-06-01',\n '2016-07-01',\n '2016-08-01',\n '2016-09-01',\n '2017-06-01',\n '2017-07-01',\n '2017-08-01',\n '2017-09-01',\n '2018-06-01',\n '2018-07-01',\n '2018-08-01',\n '2018-09-01',\n '2019-06-01',\n '2019-07-01',\n '2019-08-01',\n '2019-09-01']"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Generate NDVI image for a specific month (or year, day)"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def get_ndvi(start_date):\n start_date = ee.Date(start_date)\n end_date = start_date.advance(1, 'month')\n image = ee.ImageCollection('COPERNICUS/S2') \\\n .filterBounds(roi) \\\n .filterDate(start_date, end_date) \\\n .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) \\\n .map(maskS2clouds) \\\n .median() \\\n .clip(roi)\n \n monthly_image = image.expression(\n \" (VNIR - RED) / sqrt((VNIR + RED))\",\n {\n 'RED': image.select(\"B4\"), \n 'VNIR': image.select(\"B8\"), \n } \n )\n \n return monthly_image.rename(start_date.format('YYYY-MM-dd'))\n \n# return monthly_image.normalizedDifference(['B8', 'B4']).rename(start_date.format('YYYY-MM-dd')) ",
"execution_count": 37,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Use map() function to loop through the dates list"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "ndvi_images = ee.ImageCollection(dates.map(get_ndvi))\nndvi_bands = ndvi_images.toBands() # convert imageCollection to an image",
"execution_count": 38,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "ndvi_bands.bandNames().getInfo()",
"execution_count": 39,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 39,
"data": {
"text/plain": "['0_2016-06-01',\n '1_2016-07-01',\n '2_2016-08-01',\n '3_2016-09-01',\n '4_2017-06-01',\n '5_2017-07-01',\n '6_2017-08-01',\n '7_2017-09-01',\n '8_2018-06-01',\n '9_2018-07-01',\n '10_2018-08-01',\n '11_2018-09-01',\n '12_2019-06-01',\n '13_2019-07-01',\n '14_2019-08-01',\n '15_2019-09-01']"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "Map.addLayer(ndvi_bands, {}, 'ndvi')",
"execution_count": 40,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "image = NEDelevation.addBands(slope).addBands(aspect).addBands(ndvi_bands)",
"execution_count": 41,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# scale = 30 # spatial resolution to calculate statistics\n\n# This can take several minutes\n\n# statistics type can be: MEAN, MAXIMUM, MEDIAN, MINIMUM, STD, MIN_MAX, SUM\n\n# geemap.zonal_statistics(image, roi, 'mean.csv', statistics_type='MEAN', scale=30)",
"execution_count": 42,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "scale = 1000\nmean = image.reduceRegions(collection=roi, reducer=ee.Reducer.mean(), scale=scale)\nstd = image.reduceRegions(collection=roi, reducer=ee.Reducer.stdDev(), scale=scale)\nminimum = image.reduceRegions(collection=roi, reducer=ee.Reducer.min(), scale=scale)\nmaximum = image.reduceRegions(collection=roi, reducer=ee.Reducer.max(), scale=scale)",
"execution_count": 43,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "result = mean.merge(std).merge(minimum).merge(maximum)\nresult_info = result.getInfo()",
"execution_count": 44,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# result_info",
"execution_count": 45,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "result_info['columns']",
"execution_count": 46,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 46,
"data": {
"text/plain": "{'0_2016-06-01': 'Float',\n '10_2018-08-01': 'Float',\n '11_2018-09-01': 'Float',\n '12_2019-06-01': 'Float',\n '13_2019-07-01': 'Float',\n '14_2019-08-01': 'Float',\n '15_2019-09-01': 'Float',\n '1_2016-07-01': 'Float',\n '2_2016-08-01': 'Float',\n '3_2016-09-01': 'Float',\n '4_2017-06-01': 'Float',\n '5_2017-07-01': 'Float',\n '6_2017-08-01': 'Float',\n '7_2017-09-01': 'Float',\n '8_2018-06-01': 'Float',\n '9_2018-07-01': 'Float',\n 'ALAND': 'Long',\n 'AWATER': 'Long',\n 'DIVISION': 'String',\n 'FUNCSTAT': 'String',\n 'GEOID': 'String',\n 'INTPTLAT': 'String',\n 'INTPTLON': 'String',\n 'LSAD': 'String',\n 'MTFCC': 'String',\n 'NAME': 'String',\n 'REGION': 'String',\n 'STATEFP': 'String',\n 'STATENS': 'String',\n 'STUSPS': 'String',\n 'aspect': 'Float',\n 'elevation': 'Float',\n 'slope': 'Float',\n 'system:index': 'String'}"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import geopandas as gpd",
"execution_count": 47,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "df = gpd.GeoDataFrame.from_features(result_info['features'])",
"execution_count": 48,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "df",
"execution_count": 49,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 49,
"data": {
"text/plain": " geometry 0_2016-06-01 \\\n0 GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1... 0.284564 \n1 POLYGON ((-75.78902 39.65843, -75.78898 39.658... 0.260565 \n2 GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1... 0.194916 \n3 POLYGON ((-75.78902 39.65843, -75.78898 39.658... 0.169765 \n4 GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1... -0.136002 \n5 POLYGON ((-75.78902 39.65843, -75.78898 39.658... -0.126619 \n6 GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1... 0.621621 \n7 POLYGON ((-75.78902 39.65843, -75.78898 39.658... 0.557600 \n\n 10_2018-08-01 11_2018-09-01 12_2019-06-01 13_2019-07-01 14_2019-08-01 \\\n0 0.283944 0.219309 0.316866 0.321156 0.286520 \n1 0.286552 0.184872 0.273753 0.290885 0.278545 \n2 0.202376 0.185845 0.200007 0.199094 0.190424 \n3 0.196017 0.181372 0.184979 0.195976 0.193160 \n4 -0.211282 -0.300452 -0.138128 -0.119868 -0.126583 \n5 -0.121586 -0.452180 -0.125007 -0.100369 -0.110864 \n6 0.653156 0.623258 0.682850 0.628371 0.599217 \n7 0.617467 0.623408 0.614060 0.590652 0.612473 \n\n 15_2019-09-01 1_2016-07-01 2_2016-08-01 ... LSAD MTFCC NAME \\\n0 0.256468 0.282624 0.341025 ... 00 G4000 Maryland \n1 0.249326 0.246035 0.175100 ... 00 G4000 Delaware \n2 0.173537 0.191998 0.152416 ... 00 G4000 Maryland \n3 0.180049 0.179679 0.196984 ... 00 G4000 Delaware \n4 -0.144259 -0.164000 -0.124886 ... 00 G4000 Maryland \n5 -0.117901 -0.105460 -0.116451 ... 00 G4000 Delaware \n6 0.605813 0.644847 0.624291 ... 00 G4000 Maryland \n7 0.594031 0.602972 0.613512 ... 00 G4000 Delaware \n\n REGION STATEFP STATENS STUSPS aspect elevation slope \n0 3 24 01714934 MD 159.739159 104.544272 0.859940 \n1 3 10 01779781 DE 136.077499 11.554610 0.146803 \n2 3 24 01714934 MD 106.560749 185.224468 1.530845 \n3 3 10 01779781 DE 110.884244 15.832690 0.285621 \n4 3 24 01714934 MD 0.000000 -0.090338 0.000000 \n5 3 10 01779781 DE 0.000000 -0.042177 0.000000 \n6 3 24 01714934 MD 360.000000 950.324707 16.734945 \n7 3 10 01779781 DE 359.950195 131.447540 4.026423 \n\n[8 rows x 34 columns]",
"text/html": "<div>\n<style scoped>\n .dataframe tbody tr th:only-of-type {\n vertical-align: middle;\n }\n\n .dataframe tbody tr th {\n vertical-align: top;\n }\n\n .dataframe thead th {\n text-align: right;\n }\n</style>\n<table border=\"1\" class=\"dataframe\">\n <thead>\n <tr style=\"text-align: right;\">\n <th></th>\n <th>geometry</th>\n <th>0_2016-06-01</th>\n <th>10_2018-08-01</th>\n <th>11_2018-09-01</th>\n <th>12_2019-06-01</th>\n <th>13_2019-07-01</th>\n <th>14_2019-08-01</th>\n <th>15_2019-09-01</th>\n <th>1_2016-07-01</th>\n <th>2_2016-08-01</th>\n <th>...</th>\n <th>LSAD</th>\n <th>MTFCC</th>\n <th>NAME</th>\n <th>REGION</th>\n <th>STATEFP</th>\n <th>STATENS</th>\n <th>STUSPS</th>\n <th>aspect</th>\n <th>elevation</th>\n <th>slope</th>\n </tr>\n </thead>\n <tbody>\n <tr>\n <th>0</th>\n <td>GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1...</td>\n <td>0.284564</td>\n <td>0.283944</td>\n <td>0.219309</td>\n <td>0.316866</td>\n <td>0.321156</td>\n <td>0.286520</td>\n <td>0.256468</td>\n <td>0.282624</td>\n <td>0.341025</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Maryland</td>\n <td>3</td>\n <td>24</td>\n <td>01714934</td>\n <td>MD</td>\n <td>159.739159</td>\n <td>104.544272</td>\n <td>0.859940</td>\n </tr>\n <tr>\n <th>1</th>\n <td>POLYGON ((-75.78902 39.65843, -75.78898 39.658...</td>\n <td>0.260565</td>\n <td>0.286552</td>\n <td>0.184872</td>\n <td>0.273753</td>\n <td>0.290885</td>\n <td>0.278545</td>\n <td>0.249326</td>\n <td>0.246035</td>\n <td>0.175100</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Delaware</td>\n <td>3</td>\n <td>10</td>\n <td>01779781</td>\n <td>DE</td>\n <td>136.077499</td>\n <td>11.554610</td>\n <td>0.146803</td>\n </tr>\n <tr>\n <th>2</th>\n <td>GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1...</td>\n <td>0.194916</td>\n <td>0.202376</td>\n <td>0.185845</td>\n <td>0.200007</td>\n <td>0.199094</td>\n <td>0.190424</td>\n <td>0.173537</td>\n <td>0.191998</td>\n <td>0.152416</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Maryland</td>\n <td>3</td>\n <td>24</td>\n <td>01714934</td>\n <td>MD</td>\n <td>106.560749</td>\n <td>185.224468</td>\n <td>1.530845</td>\n </tr>\n <tr>\n <th>3</th>\n <td>POLYGON ((-75.78902 39.65843, -75.78898 39.658...</td>\n <td>0.169765</td>\n <td>0.196017</td>\n <td>0.181372</td>\n <td>0.184979</td>\n <td>0.195976</td>\n <td>0.193160</td>\n <td>0.180049</td>\n <td>0.179679</td>\n <td>0.196984</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Delaware</td>\n <td>3</td>\n <td>10</td>\n <td>01779781</td>\n <td>DE</td>\n <td>110.884244</td>\n <td>15.832690</td>\n <td>0.285621</td>\n </tr>\n <tr>\n <th>4</th>\n <td>GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1...</td>\n <td>-0.136002</td>\n <td>-0.211282</td>\n <td>-0.300452</td>\n <td>-0.138128</td>\n <td>-0.119868</td>\n <td>-0.126583</td>\n <td>-0.144259</td>\n <td>-0.164000</td>\n <td>-0.124886</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Maryland</td>\n <td>3</td>\n <td>24</td>\n <td>01714934</td>\n <td>MD</td>\n <td>0.000000</td>\n <td>-0.090338</td>\n <td>0.000000</td>\n </tr>\n <tr>\n <th>5</th>\n <td>POLYGON ((-75.78902 39.65843, -75.78898 39.658...</td>\n <td>-0.126619</td>\n <td>-0.121586</td>\n <td>-0.452180</td>\n <td>-0.125007</td>\n <td>-0.100369</td>\n <td>-0.110864</td>\n <td>-0.117901</td>\n <td>-0.105460</td>\n <td>-0.116451</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Delaware</td>\n <td>3</td>\n <td>10</td>\n <td>01779781</td>\n <td>DE</td>\n <td>0.000000</td>\n <td>-0.042177</td>\n <td>0.000000</td>\n </tr>\n <tr>\n <th>6</th>\n <td>GEOMETRYCOLLECTION (LINESTRING (-77.52722 39.1...</td>\n <td>0.621621</td>\n <td>0.653156</td>\n <td>0.623258</td>\n <td>0.682850</td>\n <td>0.628371</td>\n <td>0.599217</td>\n <td>0.605813</td>\n <td>0.644847</td>\n <td>0.624291</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Maryland</td>\n <td>3</td>\n <td>24</td>\n <td>01714934</td>\n <td>MD</td>\n <td>360.000000</td>\n <td>950.324707</td>\n <td>16.734945</td>\n </tr>\n <tr>\n <th>7</th>\n <td>POLYGON ((-75.78902 39.65843, -75.78898 39.658...</td>\n <td>0.557600</td>\n <td>0.617467</td>\n <td>0.623408</td>\n <td>0.614060</td>\n <td>0.590652</td>\n <td>0.612473</td>\n <td>0.594031</td>\n <td>0.602972</td>\n <td>0.613512</td>\n <td>...</td>\n <td>00</td>\n <td>G4000</td>\n <td>Delaware</td>\n <td>3</td>\n <td>10</td>\n <td>01779781</td>\n <td>DE</td>\n <td>359.950195</td>\n <td>131.447540</td>\n <td>4.026423</td>\n </tr>\n </tbody>\n</table>\n<p>8 rows × 34 columns</p>\n</div>"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"_draft": {
"nbviewer_url": "https://gist.github.com/365a74a0fd84332825d5409c3e07c32c"
},
"gist": {
"id": "365a74a0fd84332825d5409c3e07c32c",
"data": {
"description": "zonal statistics",
"public": true
}
},
"hide_input": false,
"kernelspec": {
"name": "python3",
"display_name": "Python 3",
"language": "python"
},
"language_info": {
"name": "python",
"version": "3.8.2",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": true,
"base_numbering": 1,
"title_cell": "Table of Contents",
"title_sidebar": "Table of Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": true
},
"varInspector": {
"window_display": false,
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"library": "var_list.py",
"delete_cmd_prefix": "del ",
"delete_cmd_postfix": "",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"library": "var_list.r",
"delete_cmd_prefix": "rm(",
"delete_cmd_postfix": ") ",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
]
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment