Skip to content

Instantly share code, notes, and snippets.

@giswqs
Created September 23, 2020 00:35
Show Gist options
  • Save giswqs/3ddb91a04c89bce5de186dd4648c5227 to your computer and use it in GitHub Desktop.
Save giswqs/3ddb91a04c89bce5de186dd4648c5227 to your computer and use it in GitHub Desktop.
zonal statistics
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import ee\nimport geemap",
"execution_count": 1,
"outputs": []
},
{
"metadata": {
"scrolled": false,
"trusted": true
},
"cell_type": "code",
"source": "Map = geemap.Map()\nMap",
"execution_count": 2,
"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": "da0f2305a5614d2c89ba1699554f48ff"
}
},
"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": 3,
"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": 4,
"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": 5,
"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": 6,
"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": 7,
"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": 8,
"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": 9,
"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": 10,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "dates.getInfo()",
"execution_count": 11,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 11,
"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": 12,
"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": 13,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "ndvi_bands.bandNames().getInfo()",
"execution_count": 14,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 14,
"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": 15,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "image = NEDelevation.addBands(slope).addBands(aspect).addBands(ndvi_bands)",
"execution_count": 16,
"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": 17,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "scale = 1000\nmean = image.reduceRegions(collection=roi, reducer=ee.Reducer.mean(), scale=scale)\nmean = mean.map(lambda f: f.set({'stat_type': 'mean'}))\nstd = image.reduceRegions(collection=roi, reducer=ee.Reducer.stdDev(), scale=scale)\nstd = std.map(lambda f: f.set({'stat_type': 'std'}))\nminimum = image.reduceRegions(collection=roi, reducer=ee.Reducer.min(), scale=scale)\nminimum = minimum.map(lambda f: f.set({'stat_type': 'min'}))\nmaximum = image.reduceRegions(collection=roi, reducer=ee.Reducer.max(), scale=scale)\nmaximum = maximum.map(lambda f: f.set({'stat_type': 'max'}))",
"execution_count": 18,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "result = mean.merge(std).merge(minimum).merge(maximum)\nresult_info = result.getInfo()",
"execution_count": 19,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# result_info",
"execution_count": 20,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "result_info['columns']",
"execution_count": 21,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 21,
"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 'stat_type': 'String',\n 'system:index': 'String'}"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import geopandas as gpd",
"execution_count": 22,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "df = gpd.GeoDataFrame.from_features(result_info['features'])",
"execution_count": 23,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "df",
"execution_count": 24,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 24,
"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 ... MTFCC NAME REGION \\\n0 0.256468 0.282624 0.341025 ... G4000 Maryland 3 \n1 0.249326 0.246035 0.175100 ... G4000 Delaware 3 \n2 0.173537 0.191998 0.152416 ... G4000 Maryland 3 \n3 0.180049 0.179679 0.196984 ... G4000 Delaware 3 \n4 -0.144259 -0.164000 -0.124886 ... G4000 Maryland 3 \n5 -0.117901 -0.105460 -0.116451 ... G4000 Delaware 3 \n6 0.605813 0.644847 0.624291 ... G4000 Maryland 3 \n7 0.594031 0.602972 0.613512 ... G4000 Delaware 3 \n\n STATEFP STATENS STUSPS aspect elevation slope stat_type \n0 24 01714934 MD 159.739159 104.544272 0.859940 mean \n1 10 01779781 DE 136.077499 11.554610 0.146803 mean \n2 24 01714934 MD 106.560749 185.224468 1.530845 std \n3 10 01779781 DE 110.884244 15.832690 0.285621 std \n4 24 01714934 MD 0.000000 -0.090338 0.000000 min \n5 10 01779781 DE 0.000000 -0.042177 0.000000 min \n6 24 01714934 MD 360.000000 950.324707 16.734945 max \n7 10 01779781 DE 359.950195 131.447540 4.026423 max \n\n[8 rows x 35 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>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 <th>stat_type</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>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 <td>mean</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>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 <td>mean</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>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 <td>std</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>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 <td>std</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>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 <td>min</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>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 <td>min</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>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 <td>max</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>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 <td>max</td>\n </tr>\n </tbody>\n</table>\n<p>8 rows × 35 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