Skip to content

Instantly share code, notes, and snippets.

@giswqs

giswqs/Untitled.ipynb

Last active September 22, 2020 21:38
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save giswqs/365a74a0fd84332825d5409c3e07c32c to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save giswqs/365a74a0fd84332825d5409c3e07c32c to your computer and use it in GitHub Desktop.
Download ZIP
zonal statistics
Display the source blob
Display the rendered blob
Raw
Untitled.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
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