Calcular las áreas de los espejos de agua a paso anual usando Google Earth Engine

Uno de los peligros por el cambio climático es el deshielo de los glaciares que permite la aparición de lagunas de manera paulatina, así como el incremento del volumen de sus aguas, y con ello el nivel de peligrosidad, ya que podrían desencadenar aluviones, perjudiciales para las poblaciones cercanas.

La temperatura de la Tierra ha aumentado en 1.1 ºC desde el periodo anterior a la Revolución Industrial y no deja de subir, al punto en que el último lustro (2015-2019) podría convertirse en el más cálido desde que hay registros (según datos de la Organización Meteorológica Mundial, publicados en diciembre del 2019).

Este incremento en la temperatura se plasma en las cordilleras del Perú, donde hay varios nevados y muchos de ellos están cerca de zonas agrícolas y pobladas que corren el peligro de acabar inundadas, donde las lagunas que se forman a su alrededor con el agua de los glaciares han incrementado su volumen. 

El uso Índice de Agua de Diferencia Normalizada Modificada o en ingles Modified Normalized Difference Water Index (MNDWI) de los satélites Landsat nos pueden ayudar a estimar el área a paso anual desde el año 1984 en adelante. El (MNDWI) utiliza bandas verdes y SWIR para mejorar las características de aguas abiertas. También disminuye las características de áreas edificadas que a menudo se correlacionan con aguas abiertas en otros índices.

En el siguiente código escrito en lenguaje Python, se hace uso de un polígono para encerrar el área de la laguna a Analizar, el código extrae las áreas año a año y lo escribe en una tabla de Excel en formato CSV. La información obtenida se puede usar para diferentes usos, como son las tendencias en la variación del volumen del cuerpo de agua o para estimar la disponibilidad hídrica de la laguna para futuros proyectos de embalses, así mismo como medida de prevención ante posibles desastres.

//Mg. Abel Carmona Arteaga

//definición de estudio

//****************definición área de estudio********************

var area_estudio = geometry;

//************//layer area de estudio***************************

var empty =ee.Image().byte();

var outline =empty.paint({

featureCollection: area_estudio,

color:1,

width:2

});

//***********************máscaras de nubes para colecciones lantsat 5, 7, 8******

var cloudMaskL457=function(image) {

 var qa=image.select('pixel_qa');

 var cloud =qa.bitwiseAnd(1 << 5)

 .and(qa.bitwiseAnd(1 << 7))

  .or(qa.bitwiseAnd(1 << 3))

  //remueve pixeles de

 var mask2 =image.mask().reduce(ee.Reducer.min());

 return image.updateMask(cloud.not()).updateMask(mask2).divide(10000).copyProperties(image, ["system:time_start"])

}

//mascara de nubes para la banda pixel_qa SR Landsat 8

function cloudMaskL8(image) {

 var cloudShadowBitMask = 1 << 3;

 var cloudsBitMask = 1 << 5;

 var qa = image.select('pixel_qa');

  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)

.and(qa.bitwiseAnd(cloudsBitMask).eq(0));

return image.updateMask(mask).divide(10000).copyProperties(image, ["system:time_start"])

}

//***********************Indices de agua*******************

function f_index_(image) {

 var mndwi=image.normalizedDifference(['green','swir1']).rename('MNDWI');

 var ndvi=image.normalizedDifference(['green','nir']).rename('MDWI');

 return image.addBands([mndwi, ndvi]).clip(area_estudio)}

   //*************** Importar colecciones Landsat************************

 var l5Bands =['B1','B2','B3','B4','B5','B6','B7','pixel_qa'];

 var l5names =['blue','green','red','nir','swir1','temp','swir2','pixel_qa'];

 var l7Bands =['B1','B2','B3','B4','B5','B6','B7','pixel_qa'];

 var l7names =['blue','green','red','nir','swir1','temp','swir2','pixel_qa'];

 var l8Bands =['B2','B3','B4','B5','B6','B10','B7','pixel_qa'];

 var l8names =['blue','green','red','nir','swir1','temp','swir2','pixel_qa'];

 //********definicion intervalo de tiempo***************

 var startyear = 1984

 var endyear = 2020;

   var startdate = ee.Date.fromYMD(startyear,1,1);

 var enddate = ee.Date.fromYMD(endyear,12,31);

   var years = ee.List.sequence(startyear, endyear);

 var months = ee.List.sequence(1,12);

   //LandSat 5 surface reflectance collection SR T1

 var ls5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR")

 .select(l5Bands, l5names)

 .filterDate(startdate,enddate)

 .map(cloudMaskL457)

 .map(f_index_)

 .filterBounds(area_estudio);

    //LandSat 7 surface reflectance collection SR T1

  var ls7 =ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")

 .select(l7Bands, l7names)

 .filterDate(startdate,enddate)

 .map(cloudMaskL457)

 .map(f_index_)

 .filterBounds(area_estudio)

    //LandSat 8 surface reflectance collection SR T1

       var ls8 =ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")

 .select(l8Bands, l8names)

 .filterDate(startdate,enddate)

 .map(cloudMaskL8)

 .map(f_index_)

 .filterBounds(area_estudio)

   //merge the dataseries

   var composites = ee.ImageCollection(ls5.merge(ls7).merge(ls8));

 var count =composites.size();

 print('numero de imagenes:',count);

   var image = ee.Image(composites.filterDate(startdate,enddate).median());

 var bandNames = image.bandNames();

 print('bands',bandNames);

   var range = composites.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])

 print('intervalo datos',ee.Date(range.get('min')), ee.Date(range.get('max')))

   //*****************Adicionando Capas**********************************

 Map.addLayer(outline, {palette:'FF0000'}, 'Borde de area de estudio');

 Map.setOptions('SATELLITE')

 Map.centerObject(area_estudio,10);

   //*************************Mascara de agua****************************

 //MNDWI maximo

 var mndwi_max = composites.select('MNDWI')

 .filterDate(startdate,enddate)

 .filter(ee.Filter.calendarRange(6,10,'month'))

 .reduce(ee.Reducer.percentile([90]))

  var mndwiMasked_max = mndwi_max.updateMask(mndwi_max.gt(0));

 var mndwiViz_max = {min:0,max:0.3068, palette:['red']};

 Map.addLayer(mndwiMasked_max.clip(area_estudio),mndwiViz_max,'MNDWI lluvioso máscara');

  //MNDWI minimo

 var mndwi_min = composites.select('MNDWI')

 .filterDate(startdate,enddate)

 .filter(ee.Filter.calendarRange(11,5,'month'))

 .reduce(ee.Reducer.percentile([90]))

  var mndwiMasked_min = mndwi_min.updateMask(mndwi_min.gt(0));

 var mndwiViz_min = {min:0,max:0.3068, palette:['cyan']};

 Map.addLayer(mndwiMasked_min.clip(area_estudio),mndwiViz_min,'MNDWI seco máscara');

  //******************Mosaico******************************************

  var imageRGB = composites.median().clip(area_estudio).visualize({bands:['red','green','blue'], min: 0.019});

   var mndwiRGB = mndwiMasked_max.clip(area_estudio).visualize({min: 0,max: 0.3068,palette:['00FFFF','#0b4a8b']

 });

   var mosaic = ee.ImageCollection([imageRGB,mndwiRGB]).mosaic();

 Map.addLayer(mosaic,{},'Mosaico',0);

   //*********************Exportando imágenes********************************

 Export.image.toDrive({

  image: mndwiMasked_max,

  description:'MNDWI_Max',

  scale: 30,

  region: area_estudio,

  folder:'Area de lagunas'

 });

   Export.image.toDrive({

  image: mndwiMasked_min,

  description:'MNDWI_Min',

  scale: 30,

  region: area_estudio,

  folder:'Area de lagunas'

 });

   //****************************cálculo de áreas****************************

 var Names =ee.List(['area de espejo de agua'])

 var totalAreas=ee.List([])

 var years =ee.List.sequence(1984,2020).getInfo()

 var month = ee.List.sequence(1,12).getInfo()

 print('Lista de años',years)

   var serie_temporal=years.map(areaCalc)

   function areaCalc(year){

  var startdate = ee.Date.fromYMD(year,1,1);

  var enddate =ee.Date.fromYMD(year,12,31).advance(4,"month");

     //Crea una nueva collecion apenas o MNDWI

     var landsat = ee.ImageCollection(composites).filterDate(startdate,enddate).select('MNDWI')

  print('Landsat' +year,landsat)

     //hacer una image

     var percMap =landsat.reduce(ee.Reducer.percentile([98])).rename('area de espejo de agua').clip(area_estudio)

  var mndwiMasked = percMap.updateMask(percMap.gt(0))

    //calculo área en Km2

    var Area =mndwiMasked

  .multiply(ee.Image.pixelArea())

  .divide(1000000)

  .reduceRegion(ee.Reducer.sum(),

  area_estudio, //suma o pixel + pixel para generar una area acumulada

  30); //suma todos los pixeles de agua

   print('area de espejo de agua' + year,Area);

      var area_total =ee.Number(Area); 

    //*****************************Creando listas array****************************

    var Areas = ee.Array(Area.get('area de espejo de agua'));

  var lista = ee.List([Areas]);

  totalAreas=totalAreas.add(Areas);

 }

 //encierra una función de cálculo de area

   //**************************Gráfico de área año a año**************************

   var chart = ui.Chart.array.values({

  array:totalAreas,

  axis: 0,

  xLabels: years,

 })

 .setSeriesNames(Names)

 .setOptions({

  title: 'Espejo de agua en Km2',

  hAxis:{title: 'serie temporal (años)', format: '####'},

  vAxis: {title: 'Area (km2)',format:'decimal'},

   // Leyenda: Nombres

  interpolateNulls: true,

  lineWidth: 1,

  pointSize: 3,

  series:{

   0:{color: 'cyan'}}//agua

  })

  .setChartType('ColumnChart');

  print(chart);