water-butoushui.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f3538c70531948aca30e947c82435415",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[30, 120], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(To…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ee\n",
    "import geemap\n",
    "#port根据自己科学工具的端口号修改\n",
    "geemap.set_proxy(port=4780)\n",
    "Map = geemap.Map(center=(30, 120), zoom=4)\n",
    "table = ee.FeatureCollection(\"users/yamiletsharon250/wuhan\")\n",
    "roi = table\n",
    "Map.addLayer(roi, {}, \"武汉\")\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#获取影像与去云\n",
    "\n",
    "#l457 去云函数\n",
    "def rmL457Cloud(image):\n",
    "    qa = image.select('pixel_qa')\n",
    "    cloud = qa.bitwiseAnd(1 << 5) \\\n",
    "                    .And(qa.bitwiseAnd(1 << 7)) \\\n",
    "                    .Or(qa.bitwiseAnd(1 << 3))\n",
    "    mask2 = image.mask().reduce(ee.Reducer.min())\n",
    "    mask3 = image.select('B1').gt(2000)\n",
    "    return image.updateMask(cloud.Not()).updateMask(mask2).updateMask(mask3.Not()) \\\n",
    "                .copyProperties(image) \\\n",
    "                .copyProperties(image, [\"system:time_start\",'system:time_end','system:footprint'])\n",
    "\n",
    "#l8 去云函数\n",
    "def rmL8Cloud(image):\n",
    "    cloudShadowBitMask = (1 << 3)\n",
    "    cloudsBitMask = (1 << 5)\n",
    "    qa = image.select('pixel_qa')\n",
    "    mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \\\n",
    "                   .And(qa.bitwiseAnd(cloudsBitMask).eq(0))\n",
    "    mask2 = image.select('B1').gt(2000)\n",
    "    return image.updateMask(mask).updateMask(mask2.Not()) \\\n",
    "                .copyProperties(image) \\\n",
    "                .copyProperties(image, [\"system:time_start\",'system:time_end'])\n",
    "                \n",
    "#定义影像执行云量筛选与去云\n",
    "l8_sr = ee.ImageCollection(\"LANDSAT/LC08/C01/T1_SR\").map(rmL8Cloud) \\\n",
    "              .filter(ee.Filter.lte('CLOUD_COVER',10))\n",
    "\n",
    "l7_sr = ee.ImageCollection(\"LANDSAT/LE07/C01/T1_SR\").map(rmL457Cloud) \\\n",
    "              .filter(ee.Filter.lte('CLOUD_COVER',10))\n",
    "\n",
    "l5_sr = ee.ImageCollection(\"LANDSAT/LT05/C01/T1_SR\").map(rmL457Cloud) \\\n",
    "              .filter(ee.Filter.lte('CLOUD_COVER',10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#计算水体方法_ARWDR，可以查阅论文，我忘了是哪篇了。。。\n",
    "#获取dem山体角,意义在于去除山体阴影\n",
    "def calcWater(image):\n",
    "    MNDWI = image.select(\"MNDWI\")\n",
    "    NDWI = image.select(\"MNDWI\")\n",
    "    NDVI = image.select(\"NDVI\")\n",
    "    EVI = image.select(\"EVI\")\n",
    "    water = NDWI.gt(0) \\\n",
    "                   .And(MNDWI.gt(0.1) \\\n",
    "                   .And(EVI.lt(0.1)) \\\n",
    "                   .And(MNDWI.gt(NDVI)) \\\n",
    "                   .Or(MNDWI.gt(EVI)))\n",
    "    dem = ee.Image(\"USGS/GMTED2010\").select(0)\n",
    "    azimuth = ee.Number(image.get('SUN_AZIMUTH'))\n",
    "    ele = ee.Number(image.get('SUN_ELEVATION'))\n",
    "    shadow = ee.Algorithms.HillShadow(dem, azimuth, ele).eq(1)\n",
    "    return image.addBands(water.rename(\"water\")).updateMask(shadow)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Landsat-5/7处理方法\n",
    "\n",
    "    #Landsat SR数据需要缩放，比例是 0.0001#\n",
    "def Landsat57scaleImage(image):\n",
    "      time_start = image.get(\"system:time_start\")\n",
    "      SUN_AZIMUTH = image.get('SOLAR_AZIMUTH_ANGLE')\n",
    "      SUN_ELEVATION = image.get('SOLAR_ZENITH_ANGLE')\n",
    "      image = image.select([\"B1\",\"B2\",\"B3\",\"B4\",\"B5\",\"B7\"])\n",
    "      image = image.divide(10000)\n",
    "      image = image.set(\"system:time_start\", time_start) \\\n",
    "                  .set(\"SUN_AZIMUTH\",SUN_AZIMUTH ) \\\n",
    "                   .set(\"SUN_ELEVATION\", SUN_ELEVATION)\n",
    "      return image\n",
    "\n",
    "    # SR数据去云#\n",
    "def Landsat57srCloudMask(image):\n",
    "      qa = image.select('pixel_qa')\n",
    "      cloudShadowBitMask = (1 << 3)\n",
    "      snowBitMask = (1 << 4)\n",
    "      cloudsBitMask = (1 << 5)\n",
    "      mask1 = qa.bitwiseAnd(cloudsBitMask).eq(0) \\\n",
    "                    .And(qa.bitwiseAnd(snowBitMask).eq(0)) \\\n",
    "                    .And(qa.bitwiseAnd(cloudShadowBitMask).eq(0))\n",
    "      mask2 = image.mask().reduce(ee.Reducer.min())\n",
    "      return image.updateMask(mask1.And(mask2))\n",
    "\n",
    "\n",
    "    #NDVI\n",
    "def Landsat57NDVI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B4\", \"B3\"]) \\\n",
    "                                   .rename(\"NDVI\"))\n",
    "    \n",
    "\n",
    "def Landsat57NDWI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B2\", \"B4\"]) \\\n",
    "                                   .rename(\"NDWI\"))\n",
    "    #MNDWI\n",
    "\n",
    "def Landsat57MNDWI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B2\", \"B5\"]) \\\n",
    "                                   .rename(\"MNDWI\"))\n",
    "\n",
    "    # EVI\n",
    "def Landsat57EVI(image):\n",
    "      evi = image.expression(\"EVI = 2.5 * (NIR - R) / (NIR + 6*R -7.5*B + 1)\", {\n",
    "        'NIR': image.select(\"B4\"),\n",
    "        'R': image.select(\"B3\"),\n",
    "        'B': image.select(\"B1\")\n",
    "      })\n",
    "      return image.addBands(evi)\n",
    "\n",
    "\n",
    "    #获取Landsat5的SR数据，并select water#\n",
    "def Landsat57getL5SRCollection(startDate, endDate, region):\n",
    "      dataset = l5_sr.filterDate(startDate, endDate) \\\n",
    "                        .filterBounds(region) \\\n",
    "                        .map(Landsat57srCloudMask) \\\n",
    "                        .map(Landsat57scaleImage) \\\n",
    "                        .map(Landsat57NDVI) \\\n",
    "                        .map(Landsat57MNDWI) \\\n",
    "                        .map(Landsat57NDWI) \\\n",
    "                        .map(Landsat57EVI) \\\n",
    "                        .map(calcWater) \\\n",
    "                        .select(\"water\")\n",
    "      return dataset\n",
    "\n",
    "    #* 获取Landsat7的SR数据#\n",
    "def Landsat57getL7SRCollection(startDate, endDate, region):\n",
    "      dataset = l7_sr.filterDate(startDate, endDate) \\\n",
    "                        .filterBounds(region) \\\n",
    "                        .map(Landsat57srCloudMask) \\\n",
    "                        .map(Landsat57scaleImage) \\\n",
    "                        .map(Landsat57NDVI) \\\n",
    "                        .map(Landsat57NDWI) \\\n",
    "                        .map(Landsat57MNDWI) \\\n",
    "                        .map(Landsat57EVI) \\\n",
    "                        .map(calcWater) \\\n",
    "                        .select(\"water\")\n",
    "      return dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Landsat 8处理方法\n",
    "\n",
    "    # Landsat8 SR数据缩放#\n",
    "def Landsat8scaleImage(image):\n",
    "      time_start = image.get(\"system:time_start\")\n",
    "      SUN_AZIMUTH = image.get('SOLAR_AZIMUTH_ANGLE')\n",
    "      SUN_ELEVATION = image.get('SOLAR_ZENITH_ANGLE')\n",
    "      image = image.select([\"B2\",\"B3\",\"B4\",\"B5\",\"B6\",\"B7\"])\n",
    "      image = image.divide(10000)\n",
    "      image = image.set(\"system:time_start\", time_start) \\\n",
    "                   .set(\"SUN_AZIMUTH\",SUN_AZIMUTH ) \\\n",
    "                   .set(\"SUN_ELEVATION\", SUN_ELEVATION)\n",
    "      return image\n",
    "\n",
    "\n",
    "    #SR数据去云#\n",
    "def Landsat8srCloudMask(image):\n",
    "      cloudShadowBitMask = (1 << 3)\n",
    "      snowBitMask = (1 << 4)\n",
    "      cloudsBitMask = (1 << 5)\n",
    "      qa = image.select('pixel_qa')\n",
    "      mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \\\n",
    "                   .And(qa.bitwiseAnd(snowBitMask).eq(0)) \\\n",
    "                   .And(qa.bitwiseAnd(cloudsBitMask).eq(0))\n",
    "      return image.updateMask(mask)\n",
    "\n",
    "\n",
    "    #NDVI\n",
    "def Landsat8NDVI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B5\", \"B4\"]) \\\n",
    "                                   .rename(\"NDVI\"))\n",
    "def Landsat8NDWI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B3\", \"B5\"]) \\\n",
    "                                   .rename(\"NDWI\"))\n",
    "\n",
    "    #MNDWI\n",
    "def Landsat8MNDWI(image):\n",
    "        return image.addBands(image.normalizedDifference([\"B3\", \"B6\"]) \\\n",
    "                                   .rename(\"MNDWI\"))\n",
    "\n",
    "\n",
    "    # EVI\n",
    "def Landsat8EVI(image):\n",
    "      evi = image.expression(\"EVI = 2.5 * (NIR - R) / (NIR + 6*R -7.5*B + 1)\", {\n",
    "        'NIR': image.select(\"B5\"),\n",
    "        'R': image.select(\"B4\"),\n",
    "        'B': image.select(\"B2\")\n",
    "      })\n",
    "      return image.addBands(evi)\n",
    "\n",
    "\n",
    "    #获取Landsat8的SR数据#\n",
    "def Landsat8getL8SRCollection(startDate, endDate, region):\n",
    "      dataset = l8_sr.filterDate(startDate, endDate) \\\n",
    "                        .filterBounds(roi) \\\n",
    "                        .map(Landsat8srCloudMask) \\\n",
    "                        .map(Landsat8scaleImage) \\\n",
    "                        .map(Landsat8NDVI) \\\n",
    "                        .map(Landsat8MNDWI) \\\n",
    "                        .map(Landsat8NDWI) \\\n",
    "                        .map(Landsat8EVI) \\\n",
    "                        .map(calcWater) \\\n",
    "                        .select(\"water\")\n",
    "      return dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2015\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f3538c70531948aca30e947c82435415",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[30, 120], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(To…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#生成逐年的水体\n",
    "def processYearWaterImage(year, region):\n",
    "    startDate = ee.Date.fromYMD(year, 1, 1)\n",
    "    endDate = ee.Date.fromYMD(year+1, 1, 1)\n",
    "    l5Water = Landsat57getL5SRCollection(startDate, endDate, region)\n",
    "    l7Water = Landsat57getL7SRCollection(startDate, endDate, region)\n",
    "    l8Water = Landsat8getL8SRCollection(startDate, endDate, region)\n",
    "    waterImgs = l5Water.merge(l7Water).merge(l8Water)\n",
    "    #计算水体的频率#\n",
    "    waterImg = waterImgs.sum() \\\n",
    "                            .divide(waterImgs.count()) \\\n",
    "                            .clip(region)\n",
    "    #提高置信度，和前面山体阴影一样，不删是因为怕报错#\n",
    "    hand = ee.ImageCollection('users/gena/global-hand/hand-100')\n",
    "    hand30 = hand.mosaic().focal_mean(0.1).rename('elevation')\n",
    "    hillShadowMask = hand30.select('elevation').lte(50)\n",
    "    dem = ee.Image(\"USGS/GMTED2010\").select(0)\n",
    "    slope = ee.Terrain.slope(dem)\n",
    "    waterImg = waterImg.updateMask(hillShadowMask).updateMask(slope.lt(20));#mask外 is NoData\n",
    "    waterImg = waterImg.unmask(0).clip(region)\n",
    "    waterImg = waterImg.where(waterImg.select([0]).gt(0).And(waterImg.select([0]).lte(0.25)),0) \\\n",
    "                      .where(waterImg.select([0]).gt(0.25).And(waterImg.select([0]).lte(0.75)),1) \\\n",
    "                      .where(waterImg.select([0]).gt(0.75).And(waterImg.select([0]).lte(1)),2) \\\n",
    "                      .toUint8() \\\n",
    "                      .clip(roi)\n",
    "    #范围从34（年份：1985）到1（年份：2018)所以是year-1985\n",
    "    butoushuiimg = ee.Image('Tsinghua/FROM-GLC/GAIA/v10').select(\"change_year_index\").eq(year-1985).clip(roi)\n",
    "    year = str(year)\n",
    "    Map.addLayer(butoushuiimg, {'min':0, 'max':0, 'opacity':1,'palette':['black']},year+\"不透水面\", False)\n",
    "    Map.addLayer(waterImg, {'min':0, 'max':2, 'opacity':1,'palette':['ffffff','red','blue']},year+\"ARWDR\", False)\n",
    "    #交集取反，留waterimg\n",
    "    difference = waterImg.updateMask(butoushuiimg.mask().Not()); \n",
    "    Map.addLayer(difference, {'min':0, 'max':2, 'opacity':1,'palette':['ffffff','red','blue']},year+\"交集\", False)\n",
    "    return waterImg\n",
    "  \n",
    "#循环导出所有的水体\n",
    "def main():\n",
    "    #开始年份和结束年份\n",
    "    imglist= []\n",
    "    year = 2015\n",
    "    endyear = 2015\n",
    "    while (year <= endyear):\n",
    "      print(year)\n",
    "      processYearWaterImage(year, roi)\n",
    "      imglist.append(processYearWaterImage(year, roi))\n",
    "      year = year + 1\n",
    "    #生成水体imagecol，方便导出\n",
    "    imglist = ee.ImageCollection(imglist)\n",
    "    #时间序列动画,有需要可以打开,不然有点卡\n",
    "    # labels = imglist.aggregate_array(\"system:index\").getInfo()\n",
    "    # Map.add_time_slider(imglist, {'min':0, 'max':2, 'opacity':0.7,'palette':['ffffff','red','blue']}, labels=labels, time_interval=3)\n",
    "\n",
    "    \n",
    "main()\n",
    "Map"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('gee')",
   "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.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5d51a8884af55115bb6d844298b2a4bf9c86da95dd33d0236e0d8eb5aff4e382"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}