From f76859504ee2ecc2776857561d218a9180f6a402 Mon Sep 17 00:00:00 2001 From: tyson-swetnam Date: Wed, 31 Oct 2018 13:10:28 -0700 Subject: [PATCH] updates .Rmd and .py --- python/gee-collection-acre-l5.py | 71 ++++ python/gee-collection-acre-l7.py | 70 ++++ python/gee-collection-acre-l8.py | 46 ++ python/gee-collection-okavango-l5.py | 71 ++++ python/gee-collection-okavango-l7.py | 70 ++++ python/gee-collection-okavango-l8.py | 45 ++ python/gee-collection-usmex-l5.py | 4 +- python/gee-collection-usmex-l7.py | 13 +- python/gee-collection-yakutia-l5.py | 71 ++++ python/gee-collection-yakutia-l7.py | 70 ++++ python/gee-collection-yakutia-l8.py | 45 ++ python/gee_l8_export.py | 45 ++ python/gee_l8_mask_export.py | 71 ++++ rmd/acre-rasters-analyses.Rmd | 422 +++++++++++++++++++ rmd/okavango-rasters-analyses.Rmd | 523 +++++++++++++++++++++++ rmd/rasters-analyses.Rmd | 499 ++++++++++++++++++++++ rmd/usmex-rasters-analyses.Rmd | 606 +++++++++++++++++++++++++++ rmd/yakutia-rasters-analyses.Rmd | 586 ++++++++++++++++++++++++++ 18 files changed, 3319 insertions(+), 9 deletions(-) create mode 100644 python/gee-collection-acre-l5.py create mode 100644 python/gee-collection-acre-l7.py create mode 100644 python/gee-collection-acre-l8.py create mode 100644 python/gee-collection-okavango-l5.py create mode 100644 python/gee-collection-okavango-l7.py create mode 100644 python/gee-collection-okavango-l8.py create mode 100644 python/gee-collection-yakutia-l5.py create mode 100644 python/gee-collection-yakutia-l7.py create mode 100644 python/gee-collection-yakutia-l8.py create mode 100644 python/gee_l8_export.py create mode 100644 python/gee_l8_mask_export.py create mode 100644 rmd/acre-rasters-analyses.Rmd create mode 100644 rmd/okavango-rasters-analyses.Rmd create mode 100644 rmd/rasters-analyses.Rmd create mode 100644 rmd/usmex-rasters-analyses.Rmd create mode 100644 rmd/yakutia-rasters-analyses.Rmd diff --git a/python/gee-collection-acre-l5.py b/python/gee-collection-acre-l5.py new file mode 100644 index 0000000..fc6412b --- /dev/null +++ b/python/gee-collection-acre-l5.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[-72.050, -9.128], + [-72.050, -9.235], + [-71.916, -9.235], + [-71.916, -9.128]]]) + +## Landsat 5 Collection +l5=ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filterDate('1998-01-01', '2011-11-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.8: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-acre-landsat5', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l5, folderName="gee-collection-acre-landsat5") diff --git a/python/gee-collection-acre-l7.py b/python/gee-collection-acre-l7.py new file mode 100644 index 0000000..db12264 --- /dev/null +++ b/python/gee-collection-acre-l7.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[-72.050, -9.128], + [-72.050, -9.235], + [-71.916, -9.235], + [-71.916, -9.128]]]) + +## Landsat 7 Collection +l7=ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filterDate('1999-01-01', '2018-09-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.7: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-acre-landsat7', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l7, folderName="gee-collection-acre-landsat7") diff --git a/python/gee-collection-acre-l8.py b/python/gee-collection-acre-l8.py new file mode 100644 index 0000000..9009b94 --- /dev/null +++ b/python/gee-collection-acre-l8.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[-72.050, -9.128], + [-72.050, -9.235], + [-71.916, -9.235], + [-71.916, -9.128]]]) + +l8=ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2013-01-01', '2018-08-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B5', 'B4']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-acre-landsat8', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(5) + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l8, folderName="gee-collection-acre-landsat8") diff --git a/python/gee-collection-okavango-l5.py b/python/gee-collection-okavango-l5.py new file mode 100644 index 0000000..5f08c8b --- /dev/null +++ b/python/gee-collection-okavango-l5.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[21.870, -18.730], + [21.870, -19.510], + [22.840, -19.510], + [22.840, -18.730]]]) + +## Landsat 5 Collection +l5=ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filterDate('1984-01-01', '1997-12-31').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.8: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-okavango-landsat5', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l5, folderName="gee-collection-okavango-landsat5") diff --git a/python/gee-collection-okavango-l7.py b/python/gee-collection-okavango-l7.py new file mode 100644 index 0000000..3a4b1aa --- /dev/null +++ b/python/gee-collection-okavango-l7.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[21.870, -18.730], + [21.870, -19.510], + [22.840, -19.510], + [22.840, -18.730]]]) + +## Landsat 7 Collection +l7=ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filterDate('1999-01-01', '2018-09-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.7: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-okavango-landsat7', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l7, folderName="gee-collection-okavango-landsat7") diff --git a/python/gee-collection-okavango-l8.py b/python/gee-collection-okavango-l8.py new file mode 100644 index 0000000..8cde8c0 --- /dev/null +++ b/python/gee-collection-okavango-l8.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +geom=ee.Geometry.Polygon( + [[[21.870, -18.730], + [21.870, -19.510], + [22.840, -19.510], + [22.840, -18.730]]]) + +l8=ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2013-01-01', '2018-08-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B5', 'B4']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-okavango-landsat8', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(5) + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l8, folderName="gee-collection-okavango-landsat8") diff --git a/python/gee-collection-usmex-l5.py b/python/gee-collection-usmex-l5.py index 0f334df..6e6b428 100644 --- a/python/gee-collection-usmex-l5.py +++ b/python/gee-collection-usmex-l5.py @@ -15,7 +15,7 @@ [-110.700, 31.429]]]) ## Landsat 5 Collection -l5=ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filterDate('1984-05-01', '2011-11-30').filterBounds(geom) +l5=ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filterDate('1998-01-01', '2011-11-30').filterBounds(geom) def exportCollectionToDrive (userCollection,folderName): userCollection2=userCollection#.map(toals) @@ -47,7 +47,7 @@ def exportImage(img): imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) ##Check for overlap - if imgp.getInfo()>=0.7: + if imgp.getInfo()>=0.8: task = ee.batch.Export.image.toDrive( image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), description = fileName, diff --git a/python/gee-collection-usmex-l7.py b/python/gee-collection-usmex-l7.py index b6ad1e7..19186f6 100644 --- a/python/gee-collection-usmex-l7.py +++ b/python/gee-collection-usmex-l7.py @@ -6,15 +6,15 @@ ee.Initialize() -## Polygon area for extraction - change AOI +## Polygon area for extraction geom=ee.Geometry.Polygon( [[[-110.700, 31.245], [-110.367, 31.245], [-110.367, 31.429], [-110.700, 31.429]]]) -## Landsat 7 Collection -l7=ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filterDate('1999-01-01', '2018-08-31').filterBounds(geom) +## Landsat 7 Collection +l7=ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filterDate('2003-05-17', '2018-09-30').filterBounds(geom) def exportCollectionToDrive (userCollection,folderName): userCollection2=userCollection#.map(toals) @@ -46,11 +46,11 @@ def exportImage(img): imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) ##Check for overlap - if imgp.getInfo()>=0.8: + if imgp.getInfo()>=0.7: task = ee.batch.Export.image.toDrive( image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), description = fileName, - folder = 'gee-collection-usmex-landsat7', + folder = 'gee-collection-usmex-landsat7-p2', maxPixels = 1e13, region = fileGeometry, scale = 30) @@ -67,5 +67,4 @@ def exportImage(img): print('Finished exporting data') print('') -## Export images to your Google Drive -exportCollectionToDrive(userCollection=l7, folderName="gee-collection-usmex-landsat7") +exportCollectionToDrive(userCollection=l7, folderName="gee-collection-usmex-landsat7-p2") diff --git a/python/gee-collection-yakutia-l5.py b/python/gee-collection-yakutia-l5.py new file mode 100644 index 0000000..3446076 --- /dev/null +++ b/python/gee-collection-yakutia-l5.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[120.080, 63.590], + [120.080, 63.160], + [121.280, 63.160], + [121.280, 63.590]]]) + +## Landsat 5 Collection +l5=ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').filterDate('1984-01-01', '1998-01-01').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.8: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-yakutia-landsat5', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l5, folderName="gee-collection-yakutia-landsat5") diff --git a/python/gee-collection-yakutia-l7.py b/python/gee-collection-yakutia-l7.py new file mode 100644 index 0000000..6682435 --- /dev/null +++ b/python/gee-collection-yakutia-l7.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[120.080, 63.590], + [120.080, 63.160], + [121.280, 63.160], + [121.280, 63.590]]]) + +## Landsat 7 Collection +l7=ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').filterDate('2003-05-17', '2018-09-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.7: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B4', 'B3']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-yakutia-landsat7-p2', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l7, folderName="gee-collection-yakutia-landsat7-p2") diff --git a/python/gee-collection-yakutia-l8.py b/python/gee-collection-yakutia-l8.py new file mode 100644 index 0000000..be6c1d5 --- /dev/null +++ b/python/gee-collection-yakutia-l8.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +geom=ee.Geometry.Polygon( + [[[120.080, 63.590], + [120.080, 63.160], + [121.280, 63.160], + [121.280, 63.590]]]) + +l8=ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2013-01-01', '2018-08-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B5', 'B4']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-yakutia-landsat8', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(5) + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l8, folderName="gee-collection-yakutia-landsat8") diff --git a/python/gee_l8_export.py b/python/gee_l8_export.py new file mode 100644 index 0000000..6d0fbfb --- /dev/null +++ b/python/gee_l8_export.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +geom=ee.Geometry.Polygon( + [[[-110.56030161232809, 31.356061922317178], + [-110.55721170754293, 31.321460260775417], + [-110.53146250099996, 31.326152757145508], + [-110.52905924172262, 31.35078452529268]]]) + +l8=ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2018-05-01', '2018-08-30').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B5', 'B4']).rename('NDVI').toFloat(), + description = fileName, + folder = 'NDVI-Test', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(5) + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l8, folderName="NDVI-Test") diff --git a/python/gee_l8_mask_export.py b/python/gee_l8_mask_export.py new file mode 100644 index 0000000..741ae25 --- /dev/null +++ b/python/gee_l8_mask_export.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import time +import oauth2client +import ee + +ee.Initialize() + +## Polygon area for extraction +## Polygon area for extraction +geom=ee.Geometry.Polygon( + [[[-110.700, 31.245], + [-110.367, 31.245], + [-110.367, 31.429], + [-110.700, 31.429]]]) + +## Landsat 8 Collection +l8=ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2013-04-11', '2018-09-11').filterBounds(geom) + +def exportCollectionToDrive (userCollection,folderName): + userCollection2=userCollection#.map(toals) + imageList = ee.List(userCollection2.toList(userCollection2.size().add(1))) + length = userCollection2.size().getInfo() + print(length) + + + def exportImage(img): + fileName = ee.String(img.get('system:index')).getInfo() + fileGeometry = geom.bounds().getInfo()['coordinates'][0] + band=ee.Image(img).bandNames().get(0) + sc=30 ##pixel size + + ## Masking + cloud=img.select('pixel_qa').bitwiseAnd(32).neq(0) + img=img.updateMask(cloud.Not()) + cloud_shadow = img.select('pixel_qa').bitwiseAnd(8).neq(0) + img=img.updateMask(cloud_shadow.Not()) + + ##Check overlap + area=ee.Feature(geom).geometry().area() + imgarea=ee.Image(img).multiply(0).add(1).clip(geom).reduceRegion( + reducer= ee.Reducer.sum().unweighted(), + geometry= geom, + maxPixels= 1e12, + tileScale= 1, + scale=30) + imgp=ee.Number(imgarea.get(band)).multiply(sc).multiply(sc).divide(area) + + ##Check for overlap + if imgp.getInfo()>=0.8: + task = ee.batch.Export.image.toDrive( + image = img.normalizedDifference(['B5', 'B4']).rename('NDVI').toFloat(), + description = fileName, + folder = 'gee-collection-usmex-landsat8', + maxPixels = 1e13, + region = fileGeometry, + scale = 30) + task.start() + + + index = 0 + while index < length: + print("Export #: " + str(index+1)) + img2export = ee.Image(imageList.get(index)) + exportImage(img2export) + index = index + 1 + time.sleep(1) # sleep for 1 seconds + + print('Finished exporting data') + print('') +exportCollectionToDrive(userCollection=l8, folderName="gee-collection-usmex-landsat8") diff --git a/rmd/acre-rasters-analyses.Rmd b/rmd/acre-rasters-analyses.Rmd new file mode 100644 index 0000000..fb599d6 --- /dev/null +++ b/rmd/acre-rasters-analyses.Rmd @@ -0,0 +1,422 @@ +--- +title: "Acre EMSI Raster Calculation" +author: "Tyson Lee Swetnam" +date: "10/22/2018" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +Instructions for installing geospatial packages are provided by [The Carpentries](https://datacarpentry.org/geospatial-workshop/setup.html) + +# Linux Dependencies +```{bash} +# sudo add-apt-repository ppa:ubuntugis +# sudo apt-get update +# sudo apt-get install libgdal-dev libgeos-dev libproj-dev +# sudo apt-get install libudunits2-dev +``` + +# R Packages +```{r message=FALSE, warning=FALSE} +# Carpentries recommended raster packages +# update.packages(checkBuilt = TRUE, ask=FALSE) +# install.packages(c('dplyr', 'ggplot2', 'raster', 'rgdal', 'rasterVis', 'remotes', 'sf')) +``` + +## A few more R Packages (just to be sure) +```{r message=FALSE, warning=FALSE} +## Install additional R Packages +# install.packages("maptools") +# install.packages("lattice") +# install.packages("RQGIS") +# install.packages("leaflet") +# install.packages("magrittr") +# install.packages("xlsx") +``` + +# Load R Libraries +```{r message=FALSE, warning=FALSE} +# Load libraries +library(rgdal) +library(grid) +library(lattice) +library(raster) +library(ggplot2) +library(RColorBrewer) +library(leaflet) +library(magrittr) +library(PerformanceAnalytics) +library(lubridate) +library(stringr) +library(reshape) +library(scales) +library(rasterVis) +library(RColorBrewer) +library(viridis) +``` + +# Import Rasters from data directories and read headers +```{r message=FALSE, warning=FALSE} +# Test for raster metadata with GDALinfo +GDALinfo("~/emsi/data/collections/acre/gee-collection-acre-landsat8/LC08_004066_20130611.tif") +``` + +## Import Raster time series for Landsats 5,7,8 +```{r message=FALSE, warning=FALSE} +# Load all rasters in acre-landsat5 +setwd("~/emsi/data/collections/acre/gee-collection-acre-landsat5/") +rlist5=list.files(getwd(), pattern=".tif$", full.names=FALSE) +for(i in rlist5) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in acre-landsat7 +setwd("~/emsi/data/collections/acre/gee-collection-acre-landsat7/") +rlist7=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist7) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in acre-landsat8 +setwd("~/emsi/data/collections/acre/gee-collection-acre-landsat8/") +rlist8=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist8) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } +``` + +```{r} +list_l5 <- ls(pattern="LT05", all.names = TRUE) +dates_l5 = as.Date(str_sub(list_l5, -8 ,-1), format="%Y%m%d") + +list_l7 <- ls(pattern="LE07", all.names = TRUE) +dates_l7 = as.Date(str_sub(list_l7, -8 ,-1), format="%Y%m%d") + +list_l8 <- ls(pattern="LC08", all.names = TRUE) +dates_l8 = as.Date(str_sub(list_l8, -8 ,-1), format="%Y%m%d") + +list_08 <- ls(pattern = "08", all.names = TRUE) +list_08 + +``` + +# August (peak greenness) EMSI calc prep +```{r} + +lall_08 <- brick(LT05_005066_19980828, + LT05_005066_19990831, + LT05_005066_20000817, + LT05_005066_20010820, + LT05_005066_20030826, + LT05_005066_20040828, + LT05_005066_20050831, + LT05_005066_20060802, + LT05_005066_20070821, + LT05_005066_20080823, + LT05_005066_20100813, + LT05_005066_20110816) + + +# LC08_005066_20130821, +# LC08_005066_20140824, +# LC08_005066_20150827, +# LC08_005066_20160829, +# LC08_005066_20170816, +# LC08_005066_20180819) + + +# Calculate mean +#l5_08_mean <- calc(l5_08, mean, na.rm=T) +#l7_08_mean <- calc(l7_08, mean, na.rm=T) +#l8_08_mean <- calc(l8_08, mean, na.rm=T) +lall_08_mean <- calc(lall_08, mean, na.rm=T) + +# Calculate sd +#l5_08_sd <- calc(l5_08, sd, na.rm=T) +#l7_08_sd <- calc(l7_08, sd, na.rm=T) +#l8_08_sd <- calc(l8_08, sd, na.rm=T) +lall_08_sd <- calc(lall_08, sd, na.rm=T) + + +#l5_08_emsi <- overlay(l5_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l7_08_emsi <- overlay(l7_08, l7_08_mean, l7_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l8_08_emsi <- overlay(l8_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +lall_08_emsi <- overlay(lall_08, lall_08_mean, lall_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +#l5_09 <- brick() + +# Landsat 8 September dates +#l8_09 <- brick(LC08_035038_20130924, +# LC08_035038_20140911, +# LC08_035038_20150930, +# LC08_035038_20160916, +# LC08_035038_20180906) +# Calculate mean +#l8_09_mean <- calc(l8_09, mean) +# Calculate sd +#l8_09_sd <- calc(l8_09, sd) +#l8_09_emsi <- overlay(l8_09, l8_09_mean, l8_09_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +# Create Leaflet Map of study area +https://rstudio.github.io/leaflet +http://leafletjs.com/ +https://www.r-bloggers.com/interactive-mapping-with-leaflet-in-r/ +https://www.color-hex.com/color-palette/19447 + +We are going to use a topo map, overlayed with a street map to show states. +To browse all the provider layers, +see http://leaflet-extras.github.io/leaflet-providers/preview/index.html + +```{r message=FALSE, warning=FALSE} +# Create custom NDVI color pallete +pal1 <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_mean), na.color = "transparent") + +pal <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_emsi[[1]]), na.color = "transparent") +``` + +```{r} +# Scene 127016 +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(lall_08_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(lall_08_mean), title = "NDVI") %>% + addRasterImage(lall_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[1]], group = "August 1998 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[2]], group = "August 1999 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[3]], group = "August 2000 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[4]], group = "August 2001 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[5]], group = "August 2003 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[6]], group = "August 2004 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[7]], group = "August 2005 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[8]], group = "August 2006 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[9]], group = "August 2007 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[10]], group = "August 2008 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[11]], group = "August 2010 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_08_emsi[[12]], group = "August 2011 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% +setView(lng = -72.050, lat = -9.128, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "August 1998 EMSI", "August 1999 EMSI", "August 2000 EMSI", "August 2001 EMSI", "August 2003 EMSI", "August 2004 EMSI", "August 2005 EMSI", "August 2006 EMSI", "August 2007 EMSI", "August 2008 EMSI", "August 2010 EMSI", "August 2011 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +#, "August 2013 EMSI", "August 2014 EMSI", "August 2015 EMSI", "August 2016 EMSI", "August 2017 EMSI", "August 2018 EMSI" +m +``` + + +```{r} +year_id <- c('LT05_005066_19980828' = "1998", + 'LT05_005066_19990831' = "1999", + 'LT05_005066_20000817' = "2000", + 'LT05_005066_20010820' = "2001", + 'LT05_005066_20030826' = "2003", + 'LT05_005066_20040828' = "2004", + 'LT05_005066_20050831' = "2005", + 'LT05_005066_20060802' = "2006", + 'LT05_005066_20070821' = "2007", + 'LT05_005066_20080823' = "2008", + 'LT05_005066_20100813' = "2010", + 'LT05_005066_20110816' = "2011") + + +year_ids <- c('layer.1'="1998", + 'layer.2'="1999", + 'layer.3'="2000", + 'layer.4'="2001", + 'layer.5'="2003", + 'layer.6'="2004", + 'layer.7'="2005", + 'layer.8'="2006", + 'layer.9'="2007", + 'layer.10'="2008", + 'layer.11'="2010", + 'layer.12'="2011") +``` + +```{r} +## Multipanel graph Summer 2000 - 2013 +lall_stack <- stack(LT05_005066_19980828, + LT05_005066_19990831, + LT05_005066_20000817, + LT05_005066_20010820, + LT05_005066_20030826, + LT05_005066_20040828, + LT05_005066_20050831, + LT05_005066_20060802, + LT05_005066_20070821, + LT05_005066_20080823, + LT05_005066_20100813, + LT05_005066_20110816) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +## Multipanel graph Augusts 2000 - 2018 +lall_stack <- stack(LT05_005066_19980828, + LT05_005066_19990831, + LT05_005066_20000817, + LT05_005066_20010820, + LT05_005066_20030826, + LT05_005066_20040828, + LT05_005066_20050831, + LT05_005066_20060802, + LT05_005066_20070821, + LT05_005066_20080823, + LT05_005066_20100813, + LT05_005066_20110816) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_acre.png', width = 22, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2010_stack_ndvi_df <- as.data.frame(lall_stack[[11]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2010_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2010_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2006_stack_ndvi_df <- as.data.frame(lall_stack[[8]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2006_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2006_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + + +```{r} +# Calculate mean +lall_stack_mean <- calc(lall_stack, mean, na.rm=T) +# Calculate sd +lall_stack_sd <- calc(lall_stack, sd, na.rm=T) +lall_stack_emsi <- overlay(lall_stack, lall_stack_mean, lall_stack_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +```{r} +lall_stack_emsi_df <- as.data.frame(lall_stack_emsi, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_timeseries_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2010_stack_emsi_df <- as.data.frame(lall_stack_emsi[[11]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2010_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2010_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2006_stack_emsi_df <- as.data.frame(lall_stack_emsi[[8]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2006_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3.2,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2006_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2016_stack_emsi_df <- as.data.frame(lall_stack_emsi[[24]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,3.2), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2016_acre.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + +```{r} +ggplot(lall_stack_emsi_df) + geom_histogram(aes(value, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + ylab("Density") + xlab("EMSI") + ggtitle("August") + + facet_wrap(~variable, labeller = as_labeller(year_ids), ncol = 6) + + theme_bw() + + theme(axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) +``` \ No newline at end of file diff --git a/rmd/okavango-rasters-analyses.Rmd b/rmd/okavango-rasters-analyses.Rmd new file mode 100644 index 0000000..257299b --- /dev/null +++ b/rmd/okavango-rasters-analyses.Rmd @@ -0,0 +1,523 @@ +--- +title: "okavango EMSI Raster Calculation" +author: "Tyson Lee Swetnam" +date: "10/13/2018" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +Instructions for installing geospatial packages are provided by [The Carpentries](https://datacarpentry.org/geospatial-workshop/setup.html) + +# Linux Dependencies +```{bash} +# sudo add-apt-repository ppa:ubuntugis +# sudo apt-get update +# sudo apt-get install libgdal-dev libgeos-dev libproj-dev +# sudo apt-get install libudunits2-dev +``` + +# R Packages +```{r message=FALSE, warning=FALSE} +# Carpentries recommended raster packages +update.packages(checkBuilt = TRUE, ask=FALSE) +install.packages(c('dplyr', 'ggplot2', 'raster', 'rgdal', 'rasterVis', 'remotes', 'sf')) +``` + +## A few more R Packages (just to be sure) +```{r message=FALSE, warning=FALSE} +## Install additional R Packages +# install.packages("maptools") +# install.packages("lattice") +# install.packages("RQGIS") +# install.packages("leaflet") +# install.packages("magrittr") +# install.packages("xlsx") +``` + +# Load R Libraries +```{r message=FALSE, warning=FALSE} +# Load libraries +library(rgdal) +library(grid) +library(lattice) +library(raster) +library(ggplot2) +library(RColorBrewer) +library(leaflet) +library(magrittr) +library(PerformanceAnalytics) +library(lubridate) +library(stringr) +library(reshape) +library(scales) +library(rasterVis) +library(RColorBrewer) +library(viridis) +``` + +# Import Rasters from data directories and read headers +```{r message=FALSE, warning=FALSE} +# Test for raster metadata with GDALinfo +GDALinfo("~/emsi/data/collections/okavango/gee-collection-okavango-landsat8/LC08_175073_20171026.tif") +``` + +## Import Raster time series for Landsats 5,7,8 +```{r message=FALSE, warning=FALSE} +# Load all rasters in okavango-landsat5 +setwd("~/emsi/data/collections/okavango/gee-collection-okavango-landsat5/") +rlist5=list.files(getwd(), pattern=".tif$", full.names=FALSE) +for(i in rlist5) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in okavango-landsat7 +setwd("~/emsi/data/collections/okavango/gee-collection-okavango-landsat7/") +rlist7=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist7) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in okavango-landsat8 +setwd("~/emsi/data/collections/okavango/gee-collection-okavango-landsat8/") +rlist8=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist8) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } +``` + +```{r} +list_l5 <- ls(pattern="LT05", all.names = TRUE) +dates_l5 = as.Date(str_sub(list_l5, -8 ,-1), format="%Y%m%d") + +list_l7 <- ls(pattern="LE07", all.names = TRUE) +dates_l7 = as.Date(str_sub(list_l7, -8 ,-1), format="%Y%m%d") + +list_l8 <- ls(pattern="LC08", all.names = TRUE) +dates_l8 = as.Date(str_sub(list_l8, -8 ,-1), format="%Y%m%d") + +list_03 <- ls(pattern = "03", all.names = TRUE) +list_03 +``` + +# March (peak greenness) EMSI calc prep +```{r} +# Landsat 5,7,8 March +lall_03 <- brick(LT05_175073_19870330, + LT05_175073_19890303, + LT05_175073_19890319, + LT05_175073_19900322, + LT05_175073_19910309, + LT05_175073_19920327, + LT05_175073_19930314, + LT05_175073_19940317, + LT05_175073_19950304, + LT05_175073_19960306, + LT05_175073_19960322, + LT05_175073_19970309, + LT05_175073_20040312, + LT05_175073_20080307, + LT05_175073_20090326, + LE07_175073_20010312, + LE07_175073_20010328, + LE07_175073_20020331, + LE07_175073_20030318, + LE07_175073_20050307, + LE07_175073_20060326, + LE07_175073_20070313, + LE07_175073_20080331, + LE07_175073_20090318, + LE07_175073_20100321, + LE07_175073_20130313, + LC08_175073_20140308, + LC08_175073_20140324, + LC08_175073_20150311, + LE07_175073_20150319, + LC08_175073_20150327, + LC08_175073_20160329, + LC08_175073_20170316, + LE07_175073_20170324, + LE07_175073_20180327, + LC08_175073_20180319) + +# LC08_175074_20140308, +# LC08_175074_20140324, +# LC08_175074_20150311, +# LC08_175074_20150327, +# LC08_175074_20160329, +# LC08_175074_20170316, +# LC08_175074_20180319, + +# LC08_127015_20140308, +# LC08_127015_20140324, +# LC08_127015_20150311, +# LC08_127015_20160313, +# LC08_127015_20170316, +# LC08_127015_20180303, +# LC08_127015_20180319, +# LC08_127016_20140308, +# LC08_127016_20140324, +# LC08_127016_20150311, +# LC08_127016_20160313, +# LC08_127016_20160329, +# LC08_127016_20180303, +# LC08_127016_20180319, +# LC08_128015_20140331, +# LC08_128015_20160304, +# LC08_128015_20170307, +# LC08_128015_20180310, +# LC08_128015_20180326, +# LC08_129015_20140306, +# LC08_129015_20140322, +# LC08_129015_20160311, +# LC08_129015_20180317, +# +# LC08_128016_20140315, +# LC08_128016_20140331, +# LC08_128016_20160304, +# LC08_128016_20180310, LC08_129016_20140306, +# LC08_129016_20160311, +# LC08_129016_20180301, +# LC08_129016_20180317, +# Calculate mean +#l5_08_mean <- calc(l5_08, mean, na.rm=T) +#l7_08_mean <- calc(l7_08, mean, na.rm=T) +#l8_08_mean <- calc(l8_08, mean, na.rm=T) +lall_03_mean <- calc(lall_03, mean, na.rm=T) + +# Calculate sd +#l5_08_sd <- calc(l5_08, sd, na.rm=T) +#l7_08_sd <- calc(l7_08, sd, na.rm=T) +#l8_08_sd <- calc(l8_08, sd, na.rm=T) +lall_03_sd <- calc(lall_03, sd, na.rm=T) + + +#l5_08_emsi <- overlay(l5_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l7_08_emsi <- overlay(l7_08, l7_08_mean, l7_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l8_08_emsi <- overlay(l8_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +lall_03_emsi <- overlay(lall_03, lall_03_mean, lall_03_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +#l5_09 <- brick() + +# Landsat 8 September dates +#l8_09 <- brick(LC08_035038_20130924, +# LC08_035038_20140911, +# LC08_035038_20150930, +# LC08_035038_20160916, +# LC08_035038_20180906) +# Calculate mean +#l8_09_mean <- calc(l8_09, mean) +# Calculate sd +#l8_09_sd <- calc(l8_09, sd) +#l8_09_emsi <- overlay(l8_09, l8_09_mean, l8_09_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +# Create Leaflet Map of study area +https://rstudio.github.io/leaflet +http://leafletjs.com/ +https://www.r-bloggers.com/interactive-mapping-with-leaflet-in-r/ +https://www.color-hex.com/color-palette/19447 + +We are going to use a topo map, overlayed with a street map to show states. +To browse all the provider layers, +see http://leaflet-extras.github.io/leaflet-providers/preview/index.html + +```{r message=FALSE, warning=FALSE} +# Create custom NDVI color pallete +pal1 <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_03_mean), na.color = "transparent") + +pal <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_03_emsi[[1]]), na.color = "transparent") +``` + +```{r} +# Okavango Early Period +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(lall_03_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(lall_03_mean), title = "NDVI") %>% + addRasterImage(lall_03_mean, group = "March Mean NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_sd, group = "March Standard Deviation NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[1]], group = "March 1987 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[2]], group = "March 1989 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[3]], group = "March 1990 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[4]], group = "March 1991 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[5]], group = "March 1992 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[6]], group = "March 1993 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[7]], group = "March 1994 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[8]], group = "March 1995 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[9]], group = "March 1996 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[10]], group = "March 1997 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[11]], group = "March 1998 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("March Mean NDVI", "March Standard Deviation NDVI", "March 1987 EMSI", "March 1989 EMSI", "March 1990 EMSI", "March 1991 EMSI", "March 1992 EMSI", "March 1993 EMSI", "March 1994 EMSI", "March 1995 EMSI", "March 1996 EMSI", "March 1997 EMSI", "March 1998 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + + +```{r} +# Okavango recent period +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(lall_03_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(lall_03_mean), title = "NDVI") %>% + addRasterImage(lall_03_mean, group = "March Mean NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_sd, group = "March Standard Deviation NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[1]], group = "March 2013 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[3]], group = "March 2014 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[5]], group = "March 2015 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_03_emsi[[6]], group = "March 2016 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% +setView(lng = 22.30, lat = -19.10, zoom = 8) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen","ESRI Satellite", "OpenStreetMap"), overlayGroups = c("March Mean NDVI", "March Standard Deviation NDVI", "March 2013 EMSI", "March 2014 EMSI", "March 2015 EMSI", "March 2016 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +year_id <- c('LT05_175073_19870330'="1987", + 'LT05_175073_19890319'="1989", + 'LT05_175073_19900322'="1990", + 'LT05_175073_19910309'="1991", + 'LT05_175073_19920327'="1992", + 'LT05_175073_19930314'="1993", + 'LT05_175073_19940317'="1994", + 'LT05_175073_19950304'="1995", + 'LT05_175073_19960322'="1996", + 'LT05_175073_19970309'="1997", + 'LE07_175073_20010312'="2001", + 'LE07_175073_20020331'="2002", + 'LE07_175073_20030318'="2003", + 'LT05_175073_20040312'="2004", + 'LE07_175073_20050307'="2005", + 'LE07_175073_20060326'="2006", + 'LE07_175073_20070313'="2007", + 'LT05_175073_20080307'="2008", + 'LT05_175073_20090326'="2009", + 'LE07_175073_20100321'="2010", + 'LE07_175073_20130313'="2013", + 'LC08_175073_20140324'="2014", + 'LC08_175073_20150311'="2015", + 'LC08_175073_20160329'="2016", + 'LC08_175073_20170316'="2017", + 'LC08_175073_20180319'="2018") +``` + +```{r} +## Multipanel graph Augusts 1987 - 2018 +lall_stack <- stack(LT05_175073_19870330, + LT05_175073_19890319, + LT05_175073_19900322, + LT05_175073_19910309, + LT05_175073_19920327, + LT05_175073_19930314, + LT05_175073_19940317, + LT05_175073_19950304, + LT05_175073_19960322, + LT05_175073_19970309, + LE07_175073_20010312, + LE07_175073_20020331, + LE07_175073_20030318, + LT05_175073_20040312, + LE07_175073_20050307, + LE07_175073_20060326, + LE07_175073_20070313, + LT05_175073_20080307, + LT05_175073_20090326, + LE07_175073_20100321, + LE07_175073_20130313, + LC08_175073_20140324, + LC08_175073_20150311, + LC08_175073_20160329, + LC08_175073_20170316, + LC08_175073_20180319) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1987_stack_ndvi_df <- as.data.frame(lall_stack[[1]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1987_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_1987_okavango.png', width = 12, height = 8, dpi = 300z, bg = "transparent") +``` +```{r} +l1996_stack_ndvi_df <- as.data.frame(lall_stack[[9]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1996_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_1996_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` +```{r} +l2016_stack_ndvi_df <- as.data.frame(lall_stack[[24]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2016_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + + +```{r} +# Calculate mean +lall_stack_mean <- calc(lall_stack, mean, na.rm=T) +# Calculate sd +lall_stack_sd <- calc(lall_stack, sd, na.rm=T) +lall_stack_emsi <- overlay(lall_stack, lall_stack_mean, lall_stack_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +# Rename variable layers back to year dates +year_ids <- c( + 'layer.1'="1987", + 'layer.2'="1989", + 'layer.3'="1990", + 'layer.4'="1991", + 'layer.5'="1992", + 'layer.6'="1993", + 'layer.7'="1994", + 'layer.8'="1995", + 'layer.9'="1996", + 'layer.10'="1997", + 'layer.11'="2001", + 'layer.12'="2002", + 'layer.13'="2003", + 'layer.14'="2004", + 'layer.15'="2005", + 'layer.16'="2006", + 'layer.17'="2007", + 'layer.18'="2008", + 'layer.19'="2009", + 'layer.20'="2010", + 'layer.21'="2013", + 'layer.22'="2014", + 'layer.23'="2015", + 'layer.24'="2016", + 'layer.25'="2017", + 'layer.26'="2018" + ) +``` + +```{r} +lall_stack_emsi_df <- as.data.frame(lall_stack_emsi, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_timeseries_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1987_stack_emsi_df <- as.data.frame(lall_stack_emsi[[1]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1987_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_1987_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1996_stack_emsi_df <- as.data.frame(lall_stack_emsi[[9]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1996_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3.2,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_1996_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2016_stack_emsi_df <- as.data.frame(lall_stack_emsi[[24]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,3.2), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2016_okavango.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + +```{r} +ggplot(lall_stack_emsi_df) + geom_histogram(aes(value, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + ylab("Density") + xlab("EMSI") + ggtitle("March") + + facet_wrap(~variable, labeller = as_labeller(year_ids), ncol = 6) + + theme_bw() + + theme(axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) +``` \ No newline at end of file diff --git a/rmd/rasters-analyses.Rmd b/rmd/rasters-analyses.Rmd new file mode 100644 index 0000000..0be41d1 --- /dev/null +++ b/rmd/rasters-analyses.Rmd @@ -0,0 +1,499 @@ +--- +title: "EMSI Raster Calculation" +author: "Tyson Lee Swetnam" +date: "9/29/2018" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + + +# Linux Dependencies + +Instructions for installing geospatial packages are provided by [The Carpentries](https://datacarpentry.org/geospatial-workshop/setup.html) + +```{bash} +# sudo add-apt-repository ppa:ubuntugis +# sudo apt-get update +# sudo apt-get install libgdal-dev libgeos-dev libproj-dev +# sudo apt-get install libudunits2-dev +``` + +# R Libraries + +## Install Packages +```{r message=FALSE, warning=FALSE} +# Carpentries recommended raster packages +# update.packages(checkBuilt = TRUE, ask=FALSE) +# install.packages(c('dplyr', 'ggplot2', 'raster', 'rgdal', 'rasterVis', 'remotes', 'sf')) +``` + +```{r message=FALSE, warning=FALSE} +## Install additional R Packages +# install.packages("maptools") +# install.packages("lattice") +# install.packages("RQGIS") +# install.packages("leaflet") +# install.packages("magrittr") +# install.packages("xlsx") +``` +## Load Libraries +```{r message=FALSE, warning=FALSE} +# Load libraries +library(rgdal) +library(grid) +library(lattice) +library(raster) +library(ggplot2) +library(RColorBrewer) +library(leaflet) +library(magrittr) +library(PerformanceAnalytics) +library(lubridate) +library(stringr) +library(reshape) +library(scales) +library(rasterVis) +library(RColorBrewer) +library(viridis) +``` + +# Import Rasters from data directory +```{r message=FALSE, warning=FALSE} +# Test for raster metadata with GDALinfo +GDALinfo("~/emsi/data/collections/usmex/gee-collection-usmex-landsat7/LE07_035038_20000912.tif") +``` + +## Import Raster time series for Landsats 5,7,8 +```{r message=FALSE, warning=FALSE} +# Load all rasters in usmex-landsat5 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat5/") +rlist5=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist5) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in usmex-landsat7 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat7/") +rlist7=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist7) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in usmex-landsat8 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat8/") +rlist8=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist8) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } +``` + +```{r} +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat5/") +rlist5 <- list.files(getwd(),full.names = TRUE,pattern = ".tif$") +l5_stack <- stack(rlist5) +l5_stack +``` + +```{r} +list_l5 <- ls(pattern="LT05", all.names = TRUE) +list_l5 +dates_l5 = as.Date(str_sub(list_l5, -8 ,-1), format="%Y%m%d") + +list_l7 <- ls(pattern="LE07", all.names = TRUE) +list_l7 +dates_l7 = as.Date(str_sub(list_l7, -8 ,-1), format="%Y%m%d") + +list_l8 <- ls(pattern="LC08", all.names = TRUE) +list_l8 +dates_l8 = as.Date(str_sub(list_l8, -8 ,-1), format="%Y%m%d") + +list_09 <- ls(pattern = "09", all.names = TRUE) +``` + +# August EMSI calc prep +```{r} +# Landsat 5,7,8 August scenes +l5_08 <- brick(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19950822, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19980830, + LT05_035038_19990817, + LT05_035038_20000803, + LT05_035038_20000819, + LT05_035038_20010822, + LT05_035038_20020809, + LT05_035038_20020825, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LT05_035038_20070823, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110802, + LT05_035038_20110818) + +l7_08 <- brick(LE07_035038_19990809, + LE07_035038_19990825, + LE07_035038_20020817, + LE07_035038_20030804, + LE07_035038_20040806, + LE07_035038_20040822, + LE07_035038_20050825, + LE07_035038_20060828, + LE07_035038_20070815, + LE07_035038_20070831, + LE07_035038_20080801, + LE07_035038_20090804, + LE07_035038_20100823, + LE07_035038_20110826, + LE07_035038_20120812, + LE07_035038_20160807, + LE07_035038_20180829) + +l8_08 <- brick( LC08_035038_20150813, + LC08_035038_20150829, + LC08_035038_20160815, + LC08_035038_20180805) + +lall_08 <- brick(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19950822, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19980830, + LT05_035038_19990817, + LT05_035038_20000803, + LT05_035038_20000819, + LT05_035038_20010822, + LT05_035038_20020809, + LT05_035038_20020825, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LT05_035038_20070823, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110802, + LT05_035038_20110818, + LE07_035038_19990809, + LE07_035038_19990825, + LE07_035038_20020817, + LE07_035038_20030804, + LE07_035038_20040806, + LE07_035038_20040822, + LE07_035038_20050825, + LE07_035038_20060828, + LE07_035038_20070815, + LE07_035038_20070831, + LE07_035038_20080801, + LE07_035038_20090804, + LE07_035038_20100823, + LE07_035038_20110826, + LE07_035038_20120812, + LE07_035038_20160807, + LE07_035038_20180829, + LC08_035038_20150813, + LC08_035038_20150829, + LC08_035038_20160815, + LC08_035038_20180805) + +# Calculate mean +l5_08_mean <- calc(l5_08, mean, na.rm=T) +l7_08_mean <- calc(l7_08, mean, na.rm=T) +l8_08_mean <- calc(l8_08, mean, na.rm=T) +lall_08_mean <- calc(lall_08, mean, na.rm=T) + +# Calculate sd +l5_08_sd <- calc(l5_08, sd, na.rm=T) +l7_08_sd <- calc(l7_08, sd, na.rm=T) +l8_08_sd <- calc(l8_08, sd, na.rm=T) +lall_08_sd <- calc(lall_08, sd, na.rm=T) + + +l5_08_emsi <- overlay(l5_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +l7_08_emsi <- overlay(l7_08, l7_08_mean, l7_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +l8_08_emsi <- overlay(l8_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +lall_08_emsi <- overlay(lall_08, lall_08_mean, lall_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +l5_09 <- brick() + +# Landsat 8 September dates +l8_09 <- brick(LC08_035038_20130924, + LC08_035038_20140911, + LC08_035038_20150930, + LC08_035038_20160916, + LC08_035038_20180906) +# Calculate mean +l8_09_mean <- calc(l8_09, mean) +# Calculate sd +l8_09_sd <- calc(l8_09, sd) +l8_09_emsi <- overlay(l8_09, l8_09_mean, l8_09_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +# Create Leaflet Map of study area +https://rstudio.github.io/leaflet +http://leafletjs.com/ +https://www.r-bloggers.com/interactive-mapping-with-leaflet-in-r/ +https://www.color-hex.com/color-palette/19447 + +We are going to use a topo map, overlayed with a street map to show states. +To browse all the provider layers, +see http://leaflet-extras.github.io/leaflet-providers/preview/index.html + + +```{r} +m <- leaflet() %>% + addTiles() %>% # Add default OpenStreetMap map tiles + addMarkers(lng=174.768, lat=-36.852, + popup="The birthplace of R") +m # Print the map +``` + +```{r message=FALSE, warning=FALSE} +# Create custom NDVI color pallete +pal1 <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_mean), na.color = "transparent") + +pal <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_emsi[[1]]), na.color = "transparent") +``` + +```{r} +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(l8_08_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(l8_08_mean), title = "NDVI") %>% + addRasterImage(l8_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l8_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[1]], group = "early August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[2]], group = "August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[3]], group = "August 2016 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[4]], group = "August 2018 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "early August 2015 EMSI", "August 2015 EMSI", "August 2016 EMSI", "August 2018 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +# Landsat 5 +m <- leaflet() %>% + addTiles() %>% + addLegend(pal = pal, values = values(l5_08_emsi[[1]]), title = "EMSI") %>% + addLegend(pal = pal1, values = values(l5_08_mean), title = "NDVI") %>% + addRasterImage(l5_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l5_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[1]], group = "August 1985 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[2]], group = "August 1987 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[3]], group = "August 1988 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[4]], group = "August 1989 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[5]], group = "August 1990 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[6]], group = "August 1991 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[7]], group = "August 1992 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[8]], group = "August 1993 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[9]], group = "August 1994 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "August 1985 EMSI", "August 1987 EMSI", "August 1988 EMSI", "August 1989 EMSI", "August 1990 EMSI", "August 1991 EMSI", "August 1992 EMSI", "August 1993 EMSI", "August 1994 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} + +library(raster) +library(rasterVis) +``` + +```{r} +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(l8_08_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(l8_08_mean), title = "NDVI") %>% + addRasterImage(lall_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(lall_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[1]], group = "early August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[2]], group = "August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[3]], group = "August 2016 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[4]], group = "August 2018 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "early August 2015 EMSI", "August 2015 EMSI", "August 2016 EMSI", "August 2018 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +year_id <- c('LT05_035038_19850826'="1985", + 'LT05_035038_19870816'="1987", + 'LT05_035038_19880818'="1985", + 'LT05_035038_19890821'="1989", + 'LT05_035038_19900824'="1990", + 'LT05_035038_19910811'="1991", + 'LT05_035038_19920813'="1992", + 'LT05_035038_19930816'="1993", + 'LT05_035038_19940803'="1994", + 'LT05_035038_19950806'="1995", + 'LT05_035038_19960824'="1996", + 'LT05_035038_19970827'="1997", + 'LT05_035038_19980814'="1998", + 'LT05_035038_19990817'="1999", + 'LT05_035038_20000819'="2000", + 'LT05_035038_20010822'="2001", + 'LE07_035038_20020817'="2002", + 'LT05_035038_20030812'="2003", + 'LT05_035038_20040830'="2004", + 'LT05_035038_20050801'="2005", + 'LE07_035038_20060828'="2006", + 'LT05_035038_20070823'="2007", + 'LE07_035038_20080801'="2008", + 'LT05_035038_20090828'="2009", + 'LT05_035038_20100815'="2010", + 'LT05_035038_20110818'="2011", + 'LE07_035038_20120812'="2012", + 'LC08_035038_20150813'="2015", + 'LC08_035038_20160815'="2016", + 'LC08_035038_20180805'="2018" + ) +``` + +```{r} +## Multipanel graph Augusts 1985 - 2018 +lall_stack <- stack(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19990817, + LT05_035038_20000819, + LT05_035038_20010822, + LE07_035038_20020817, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LE07_035038_20060828, + LT05_035038_20070823, + LE07_035038_20080801, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110818, + LE07_035038_20120812, + LC08_035038_20150813, + LC08_035038_20160815, + LC08_035038_20180805) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +# Calculate mean +lall_stack_mean <- calc(lall_stack, mean, na.rm=T) +# Calculate sd +lall_stack_sd <- calc(lall_stack, sd, na.rm=T) +lall_stack_emsi <- overlay(lall_stack, lall_stack_mean, lall_stack_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +# Rename variable layers back to year dates +year_ids <- c( + 'layer.1'="1985", + 'layer.2'="1987", + 'layer.3'="1988", + 'layer.4'="1989", + 'layer.5'="1990", + 'layer.6'="1991", + 'layer.7'="1992", + 'layer.8'="1993", + 'layer.9'="1994", + 'layer.10'="1995", + 'layer.11'="1996", + 'layer.12'="1997", + 'layer.13'="1998", + 'layer.14'="1999", + 'layer.15'="2000", + 'layer.16'="2001", + 'layer.17'="2002", + 'layer.18'="2003", + 'layer.19'="2004", + 'layer.20'="2005", + 'layer.21'="2006", + 'layer.22'="2007", + 'layer.23'="2008", + 'layer.24'="2009", + 'layer.25'="2010", + 'layer.26'="2011", + 'layer.27'="2012", + 'layer.28'="2015", + 'layer.29'="2016", + 'layer.30'="2018" + ) +``` + +```{r} +lall_stack_emsi_df <- as.data.frame(lall_stack_emsi, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_timeseries_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` \ No newline at end of file diff --git a/rmd/usmex-rasters-analyses.Rmd b/rmd/usmex-rasters-analyses.Rmd new file mode 100644 index 0000000..60e860d --- /dev/null +++ b/rmd/usmex-rasters-analyses.Rmd @@ -0,0 +1,606 @@ +--- +title: "US Mexico EMSI Raster Calculation" +author: "Tyson Lee Swetnam" +date: "10/09/2018" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + + +# Linux Dependencies + +Instructions for installing geospatial packages are provided by [The Carpentries](https://datacarpentry.org/geospatial-workshop/setup.html) + +```{bash} +# sudo add-apt-repository ppa:ubuntugis +# sudo apt-get update +# sudo apt-get install libgdal-dev libgeos-dev libproj-dev +# sudo apt-get install libudunits2-dev +``` + +# R Libraries + +## Install Packages +```{r message=FALSE, warning=FALSE} +# Carpentries recommended raster packages +# update.packages(checkBuilt = TRUE, ask=FALSE) +# install.packages(c('dplyr', 'ggplot2', 'raster', 'rgdal', 'rasterVis', 'remotes', 'sf')) +``` + +```{r message=FALSE, warning=FALSE} +## Install additional R Packages +# install.packages("maptools") +# install.packages("lattice") +# install.packages("RQGIS") +# install.packages("leaflet") +# install.packages("magrittr") +# install.packages("xlsx") +``` +## Load Libraries +```{r message=FALSE, warning=FALSE} +# Load libraries +library(rgdal) +library(grid) +library(lattice) +library(raster) +library(ggplot2) +library(RColorBrewer) +library(leaflet) +library(magrittr) +library(PerformanceAnalytics) +library(lubridate) +library(stringr) +library(reshape) +library(scales) +library(rasterVis) +library(RColorBrewer) +library(viridis) +``` + +# Import Rasters from data directory +```{r message=FALSE, warning=FALSE} +# Test for raster metadata with GDALinfo +GDALinfo("~/emsi/data/collections/usmex/gee-collection-usmex-landsat7/LE07_035038_20000912.tif") +``` + +## Import Raster time series for Landsats 5,7,8 +```{r message=FALSE, warning=FALSE} +# Load all rasters in usmex-landsat5 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat5/") +rlist5=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist5) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in usmex-landsat7 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat7/") +rlist7=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist7) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in usmex-landsat8 +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat8/") +rlist8=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist8) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } +``` + +```{r} +setwd("~/emsi/data/collections/usmex/gee-collection-usmex-landsat5/") +rlist5 <- list.files(getwd(),full.names = TRUE,pattern = ".tif$") +l5_stack <- stack(rlist5) +l5_stack +``` + +```{r} +list_l5 <- ls(pattern="LT05", all.names = TRUE) +list_l5 +dates_l5 = as.Date(str_sub(list_l5, -8 ,-1), format="%Y%m%d") + +list_l7 <- ls(pattern="LE07", all.names = TRUE) +list_l7 +dates_l7 = as.Date(str_sub(list_l7, -8 ,-1), format="%Y%m%d") + +list_l8 <- ls(pattern="LC08", all.names = TRUE) +list_l8 +dates_l8 = as.Date(str_sub(list_l8, -8 ,-1), format="%Y%m%d") + +list_09 <- ls(pattern = "09", all.names = TRUE) +``` + +# August EMSI calc prep +```{r} +# Landsat 5,7,8 August scenes +l5_08 <- brick(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19950822, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19980830, + LT05_035038_19990817, + LT05_035038_20000803, + LT05_035038_20000819, + LT05_035038_20010822, + LT05_035038_20020809, + LT05_035038_20020825, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LT05_035038_20070823, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110802, + LT05_035038_20110818) + +l7_08 <- brick(LE07_035038_19990809, + LE07_035038_19990825, + LE07_035038_20020817, + LE07_035038_20030804, + LE07_035038_20040806, + LE07_035038_20040822, + LE07_035038_20050825, + LE07_035038_20060828, + LE07_035038_20070815, + LE07_035038_20070831, + LE07_035038_20080801, + LE07_035038_20090804, + LE07_035038_20100823, + LE07_035038_20110826, + LE07_035038_20120812, + LE07_035038_20160807, + LE07_035038_20180829) + +l8_08 <- brick( LC08_035038_20150813, + LC08_035038_20150829, + LC08_035038_20160815, + LC08_035038_20180805) + +lall_08 <- brick(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19950822, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19980830, + LT05_035038_19990817, + LT05_035038_20000803, + LT05_035038_20000819, + LT05_035038_20010822, + LT05_035038_20020809, + LT05_035038_20020825, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LT05_035038_20070823, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110802, + LT05_035038_20110818, + LE07_035038_19990809, + LE07_035038_19990825, + LE07_035038_20020817, + LE07_035038_20030804, + LE07_035038_20040806, + LE07_035038_20040822, + LE07_035038_20050825, + LE07_035038_20060828, + LE07_035038_20070815, + LE07_035038_20070831, + LE07_035038_20080801, + LE07_035038_20090804, + LE07_035038_20100823, + LE07_035038_20110826, + LE07_035038_20120812, + LE07_035038_20160807, + LE07_035038_20180829, + LC08_035038_20150813, + LC08_035038_20150829, + LC08_035038_20160815, + LC08_035038_20180805) + +# Calculate mean +l5_08_mean <- calc(l5_08, mean, na.rm=T) +l7_08_mean <- calc(l7_08, mean, na.rm=T) +l8_08_mean <- calc(l8_08, mean, na.rm=T) +lall_08_mean <- calc(lall_08, mean, na.rm=T) + +# Calculate sd +l5_08_sd <- calc(l5_08, sd, na.rm=T) +l7_08_sd <- calc(l7_08, sd, na.rm=T) +l8_08_sd <- calc(l8_08, sd, na.rm=T) +lall_08_sd <- calc(lall_08, sd, na.rm=T) + + +l5_08_emsi <- overlay(l5_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +l7_08_emsi <- overlay(l7_08, l7_08_mean, l7_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +l8_08_emsi <- overlay(l8_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +lall_08_emsi <- overlay(lall_08, lall_08_mean, lall_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +l5_09 <- brick() + +# Landsat 8 September dates +l8_09 <- brick(LC08_035038_20130924, + LC08_035038_20140911, + LC08_035038_20150930, + LC08_035038_20160916, + LC08_035038_20180906) +# Calculate mean +l8_09_mean <- calc(l8_09, mean) +# Calculate sd +l8_09_sd <- calc(l8_09, sd) +l8_09_emsi <- overlay(l8_09, l8_09_mean, l8_09_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +# Create Leaflet Map of study area +https://rstudio.github.io/leaflet +http://leafletjs.com/ +https://www.r-bloggers.com/interactive-mapping-with-leaflet-in-r/ +https://www.color-hex.com/color-palette/19447 + +We are going to use a topo map, overlayed with a street map to show states. +To browse all the provider layers, +see http://leaflet-extras.github.io/leaflet-providers/preview/index.html + + +```{r} +m <- leaflet() %>% + addTiles() %>% # Add default OpenStreetMap map tiles + addMarkers(lng=174.768, lat=-36.852, + popup="The birthplace of R") +m # Print the map +``` + +```{r message=FALSE, warning=FALSE} +# Create custom NDVI color pallete +pal1 <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_mean), na.color = "transparent") + +pal <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_08_emsi[[1]]), na.color = "transparent") +``` + +```{r} +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(l8_08_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(l8_08_mean), title = "NDVI") %>% + addRasterImage(l8_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l8_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[1]], group = "early August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[2]], group = "August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[3]], group = "August 2016 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l8_08_emsi[[4]], group = "August 2018 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "early August 2015 EMSI", "August 2015 EMSI", "August 2016 EMSI", "August 2018 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +# Landsat 5 +m <- leaflet() %>% + addTiles() %>% + addLegend(pal = pal, values = values(l5_08_emsi[[1]]), title = "EMSI") %>% + addLegend(pal = pal1, values = values(l5_08_mean), title = "NDVI") %>% + addRasterImage(l5_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l5_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[1]], group = "August 1985 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[2]], group = "August 1987 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[3]], group = "August 1988 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[4]], group = "August 1989 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[5]], group = "August 1990 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[6]], group = "August 1991 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[7]], group = "August 1992 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[8]], group = "August 1993 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(l5_08_emsi[[9]], group = "August 1994 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "August 1985 EMSI", "August 1987 EMSI", "August 1988 EMSI", "August 1989 EMSI", "August 1990 EMSI", "August 1991 EMSI", "August 1992 EMSI", "August 1993 EMSI", "August 1994 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} + +library(raster) +library(rasterVis) +``` + +```{r} +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(l8_08_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(l8_08_mean), title = "NDVI") %>% + addRasterImage(lall_08_mean, group = "August Mean NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(lall_08_sd, group = "August Standard Deviation NDVI", colors = pal1, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[1]], group = "early August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[2]], group = "August 2015 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[3]], group = "August 2016 EMSI", colors = pal, opacity = 1.0) %>% + addRasterImage(lall_08_emsi[[4]], group = "August 2018 EMSI", colors = pal, opacity = 1.0) %>% +setView(lng = -110.55, lat = 31.3, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("August Mean NDVI", "August Standard Deviation NDVI", "early August 2015 EMSI", "August 2015 EMSI", "August 2016 EMSI", "August 2018 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +year_id <- c('LT05_035038_19850826'="1985", + 'LT05_035038_19870816'="1987", + 'LT05_035038_19880818'="1985", + 'LT05_035038_19890821'="1989", + 'LT05_035038_19900824'="1990", + 'LT05_035038_19910811'="1991", + 'LT05_035038_19920813'="1992", + 'LT05_035038_19930816'="1993", + 'LT05_035038_19940803'="1994", + 'LT05_035038_19950806'="1995", + 'LT05_035038_19960824'="1996", + 'LT05_035038_19970827'="1997", + 'LT05_035038_19980814'="1998", + 'LT05_035038_19990817'="1999", + 'LT05_035038_20000819'="2000", + 'LT05_035038_20010822'="2001", + 'LE07_035038_20020817'="2002", + 'LT05_035038_20030812'="2003", + 'LT05_035038_20040830'="2004", + 'LT05_035038_20050801'="2005", + 'LE07_035038_20060828'="2006", + 'LT05_035038_20070823'="2007", + 'LE07_035038_20080801'="2008", + 'LT05_035038_20090828'="2009", + 'LT05_035038_20100815'="2010", + 'LT05_035038_20110818'="2011", + 'LE07_035038_20120812'="2012", + 'LC08_035038_20150813'="2015", + 'LC08_035038_20160815'="2016", + 'LC08_035038_20180805'="2018" + ) +``` + +```{r} +## Multipanel graph Augusts 1985 - 2018 +lall_stack <- stack(LT05_035038_19850826, + LT05_035038_19870816, + LT05_035038_19880818, + LT05_035038_19890821, + LT05_035038_19900824, + LT05_035038_19910811, + LT05_035038_19920813, + LT05_035038_19930816, + LT05_035038_19940803, + LT05_035038_19950806, + LT05_035038_19960824, + LT05_035038_19970827, + LT05_035038_19980814, + LT05_035038_19990817, + LT05_035038_20000819, + LT05_035038_20010822, + LE07_035038_20020817, + LT05_035038_20030812, + LT05_035038_20040830, + LT05_035038_20050801, + LE07_035038_20060828, + LT05_035038_20070823, + LE07_035038_20080801, + LT05_035038_20090828, + LT05_035038_20100815, + LT05_035038_20110818, + LE07_035038_20120812, + LC08_035038_20150813, + LC08_035038_20160815, + LC08_035038_20180805) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1994_stack_ndvi_df <- as.data.frame(lall_stack[[9]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1994_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_1994_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` +```{r} +l2001_stack_ndvi_df <- as.data.frame(lall_stack[[16]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2001_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2001_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` +```{r} +l2016_stack_ndvi_df <- as.data.frame(lall_stack[[29]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2016_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + + +```{r} +# Calculate mean +lall_stack_mean <- calc(lall_stack, mean, na.rm=T) +# Calculate sd +lall_stack_sd <- calc(lall_stack, sd, na.rm=T) +lall_stack_emsi <- overlay(lall_stack, lall_stack_mean, lall_stack_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +# Rename variable layers back to year dates +year_ids <- c( + 'layer.1'="1985", + 'layer.2'="1987", + 'layer.3'="1988", + 'layer.4'="1989", + 'layer.5'="1990", + 'layer.6'="1991", + 'layer.7'="1992", + 'layer.8'="1993", + 'layer.9'="1994", + 'layer.10'="1995", + 'layer.11'="1996", + 'layer.12'="1997", + 'layer.13'="1998", + 'layer.14'="1999", + 'layer.15'="2000", + 'layer.16'="2001", + 'layer.17'="2002", + 'layer.18'="2003", + 'layer.19'="2004", + 'layer.20'="2005", + 'layer.21'="2006", + 'layer.22'="2007", + 'layer.23'="2008", + 'layer.24'="2009", + 'layer.25'="2010", + 'layer.26'="2011", + 'layer.27'="2012", + 'layer.28'="2015", + 'layer.29'="2016", + 'layer.30'="2018" + ) +``` + +```{r} +lall_stack_emsi_df <- as.data.frame(lall_stack_emsi, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_timeseries_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1994_stack_emsi_df <- as.data.frame(lall_stack_emsi[[9]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1994_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_1994_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2001_stack_emsi_df <- as.data.frame(lall_stack_emsi[[16]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2001_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2001_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` +```{r} +l2016_stack_emsi_df <- as.data.frame(lall_stack_emsi[[29]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2016_grassmx.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + +```{r} +ggplot(lall_stack_emsi_df) + geom_histogram(aes(value, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + ylab("Density") + xlab("EMSI") + ggtitle("August") + + facet_wrap(~variable, labeller = as_labeller(year_ids), ncol = 6) + + theme_bw() + + theme(axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) +``` \ No newline at end of file diff --git a/rmd/yakutia-rasters-analyses.Rmd b/rmd/yakutia-rasters-analyses.Rmd new file mode 100644 index 0000000..28e9824 --- /dev/null +++ b/rmd/yakutia-rasters-analyses.Rmd @@ -0,0 +1,586 @@ +--- +title: "Yakutia EMSI Raster Calculation" +author: "Tyson Lee Swetnam" +date: "10/13/2018" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +Instructions for installing geospatial packages are provided by [The Carpentries](https://datacarpentry.org/geospatial-workshop/setup.html) + +# Linux Dependencies +```{bash} +# sudo add-apt-repository ppa:ubuntugis +# sudo apt-get update +# sudo apt-get install libgdal-dev libgeos-dev libproj-dev +# sudo apt-get install libudunits2-dev +``` + +# R Packages +```{r message=FALSE, warning=FALSE} +# Carpentries recommended raster packages +# update.packages(checkBuilt = TRUE, ask=FALSE) +# install.packages(c('dplyr', 'ggplot2', 'raster', 'rgdal', 'rasterVis', 'remotes', 'sf')) +``` + +## A few more R Packages (just to be sure) +```{r message=FALSE, warning=FALSE} +## Install additional R Packages +# install.packages("maptools") +# install.packages("lattice") +# install.packages("RQGIS") +# install.packages("leaflet") +# install.packages("magrittr") +# install.packages("xlsx") +``` + +# Load R Libraries +```{r message=FALSE, warning=FALSE} +# Load libraries +library(rgdal) +library(grid) +library(lattice) +library(raster) +library(ggplot2) +library(RColorBrewer) +library(leaflet) +library(magrittr) +library(PerformanceAnalytics) +library(lubridate) +library(stringr) +library(reshape) +library(scales) +library(rasterVis) +library(RColorBrewer) +library(viridis) +``` + +# Import Rasters from data directories and read headers +```{r message=FALSE, warning=FALSE} +# Test for raster metadata with GDALinfo +GDALinfo("~/emsi/data/collections/yakutia/gee-collection-yakutia-landsat8/LC08_175073_20171026.tif") +``` + +## Import Raster time series for Landsats 5,7,8 +```{r message=FALSE, warning=FALSE} +# Load all rasters in yakutia-landsat5 +setwd("~/emsi/data/collections/yakutia/gee-collection-yakutia-landsat5/") +rlist5=list.files(getwd(), pattern=".tif$", full.names=FALSE) +for(i in rlist5) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in yakutia-landsat7 +setwd("~/emsi/data/collections/yakutia/gee-collection-yakutia-landsat7/") +rlist7=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist7) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } + +# Load all rasters in yakutia-landsat8 +setwd("~/emsi/data/collections/yakutia/gee-collection-yakutia-landsat8/") +rlist8=list.files(getwd(), pattern="tif$", full.names=FALSE) +for(i in rlist8) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) } +``` + +```{r} +list_l5 <- ls(pattern="LT05", all.names = TRUE) +dates_l5 = as.Date(str_sub(list_l5, -8 ,-1), format="%Y%m%d") + +list_l7 <- ls(pattern="LE07", all.names = TRUE) +dates_l7 = as.Date(str_sub(list_l7, -8 ,-1), format="%Y%m%d") + +list_l8 <- ls(pattern="LC08", all.names = TRUE) +dates_l8 = as.Date(str_sub(list_l8, -8 ,-1), format="%Y%m%d") + +list_06 <- ls(pattern = "06", all.names = TRUE) +list_06 + +list_07 <- ls(pattern = "07", all.names = TRUE) +list_07 +``` + +# June (peak greenness) EMSI calc prep +```{r} + +lall_06 <- brick(LT05_128016_19990728, + LT05_128016_20010701, + LT05_128016_20010717, + LT05_128016_20020704, + LT05_128016_20020720, + LE07_128016_20030715, + LT05_128016_20040725, + LE07_128016_20050704, + LE07_128016_20060520, + LE07_128016_20060723, + LT05_128016_20060731, + LT05_128016_20070515, + LE07_128016_20070726, + LE07_128016_20080626, + LE07_128016_20080813, + LE07_128016_20090715, + LE07_128016_20090917, + LE07_128016_20100819, + LE07_128016_20110603, + LT05_128016_20110713, + LE07_128016_20120707, + LE07_128016_20120723, + LE07_128016_20130726, + LE07_128016_20140611, + LC08_128016_20140619, + LC08_128016_20150606, + LC08_128016_20150622, + LE07_128016_20160803, + LE07_128016_20160920, + LC08_128016_20160608, + LC08_128016_20160624, + LC08_128016_20170611, + LC08_128016_20170627, + LE07_128016_20170705, + LE07_128016_20170721, + LC08_128016_20180614, + LC08_128016_20180630, + LE07_128016_20180521, + LE07_128016_20180606 + ) + + +# Landsat 5,7,8 June +lall_06 <- brick(LT05_127016_20000605, + LT05_127016_20010624, + LT05_127016_20020627, + LE07_127016_20050611, + LT05_127016_20060910, + LT05_127016_20070812, + LT05_127016_20090630, + LT05_127016_20100601, + LE07_127016_20120630, + LC08_127016_20130609) + + +# Calculate mean +#l5_08_mean <- calc(l5_08, mean, na.rm=T) +#l7_08_mean <- calc(l7_08, mean, na.rm=T) +#l8_08_mean <- calc(l8_08, mean, na.rm=T) +lall_06_mean <- calc(lall_06, mean, na.rm=T) + +# Calculate sd +#l5_08_sd <- calc(l5_08, sd, na.rm=T) +#l7_08_sd <- calc(l7_08, sd, na.rm=T) +#l8_08_sd <- calc(l8_08, sd, na.rm=T) +lall_06_sd <- calc(lall_06, sd, na.rm=T) + + +#l5_08_emsi <- overlay(l5_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l7_08_emsi <- overlay(l7_08, l7_08_mean, l7_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +#l8_08_emsi <- overlay(l8_08, l5_08_mean, l5_08_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) + +lall_06_emsi <- overlay(lall_06, lall_06_mean, lall_06_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +#l5_09 <- brick() + +# Landsat 8 September dates +#l8_09 <- brick(LC08_035038_20130924, +# LC08_035038_20140911, +# LC08_035038_20150930, +# LC08_035038_20160916, +# LC08_035038_20180906) +# Calculate mean +#l8_09_mean <- calc(l8_09, mean) +# Calculate sd +#l8_09_sd <- calc(l8_09, sd) +#l8_09_emsi <- overlay(l8_09, l8_09_mean, l8_09_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + + +# Create Leaflet Map of study area +https://rstudio.github.io/leaflet +http://leafletjs.com/ +https://www.r-bloggers.com/interactive-mapping-with-leaflet-in-r/ +https://www.color-hex.com/color-palette/19447 + +We are going to use a topo map, overlayed with a street map to show states. +To browse all the provider layers, +see http://leaflet-extras.github.io/leaflet-providers/preview/index.html + +```{r message=FALSE, warning=FALSE} +# Create custom NDVI color pallete +pal1 <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_06_mean), na.color = "transparent") + +pal <- colorNumeric(c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), values(lall_06_emsi[[1]]), na.color = "transparent") +``` + +```{r} +# Scene 127016 +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(lall_06_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(lall_06_mean), title = "NDVI") %>% + addRasterImage(lall_06_mean, group = "June Mean NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_sd, group = "June Standard Deviation NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[1]], group = "June 2000 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[2]], group = "June 2001 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[3]], group = "June 2002 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[4]], group = "June 2005 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[5]], group = "July 2006 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[6]], group = "July 2007 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[7]], group = "July 2009 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[8]], group = "July 2010 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[9]], group = "July 2012 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[10]], group = "July 2013 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% +setView(lng = 120.3300, lat = 63.4300, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("June Mean NDVI", "June Standard Deviation NDVI", "June 2000 EMSI", "June 2001 EMSI", "June 2002 EMSI", "June 2005 EMSI", "June 2006 EMSI", "July 2007 EMSI", "July 2009 EMSI", "July 2010 EMSI", "July 2012 EMSI", "July 2013 EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + +```{r} +# Scene 127016 +m <- leaflet() %>% + addTiles() %>% + #addLegend(pal = pal, values = values(lall_06_emsi[[1]]), title = "EMSI") %>% + #addLegend(pal = pal1, values = values(lall_06_mean), title = "NDVI") %>% + addRasterImage(lall_06_mean, group = "June Mean NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_sd, group = "June Standard Deviation NDVI", colors = pal1, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[21]], group = "June 2013 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[22]], group = "June 2014a EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[23]], group = "June 2014b EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[24]], group = "June 2015 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[25]], group = "July 2015a EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[26]], group = "July 2015b EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[27]], group = "July 2016 EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[28]], group = "July 2017a EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% + addRasterImage(lall_06_emsi[[29]], group = "July 2017b EMSI", colors = pal, opacity = 1.0, maxBytes = 32 * 1024 * 1024) %>% +setView(lng = 120.3300, lat = 63.4300, zoom = 12) %>% +addProviderTiles("Stamen.Toner", group = "Stamen") %>% +addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite", options = providerTileOptions(opacity = 0.66, transparent = TRUE)) %>% +addProviderTiles("OpenStreetMap.Mapnik", group = "OpenStreetMap") %>% +#layers control panel +addLayersControl(baseGroups = c("Stamen", "ESRI Satellite", "OpenStreetMap"), overlayGroups = c("June Mean NDVI", "June Standard Deviation NDVI", "June 2013 EMSI", "June 2014a EMSI", "June 2014b EMSI", "June 2015 EMSI", "July 2015a EMSI", "July 2015b EMSI", "July 2016 EMSI", "July 2017a EMSI", "July 2017b EMSI"), options = layersControlOptions(collapsed = TRUE)) + +m +``` + + +```{r} + +year_id <- c('LT05_128016_19990728' = "1999", + 'LT05_128016_20010701' = "2001", + 'LT05_128016_20010717' = "2001", + 'LT05_128016_20020704' = "2002", + 'LT05_128016_20020720' = "2002", + 'LE07_128016_20030715' = "2003", + 'LT05_128016_20040725' = "2004", + 'LE07_128016_20050704' = "2005", + 'LE07_128016_20060520' = "2006", + 'LE07_128016_20060723' = "2006", + 'LT05_128016_20060731' = "2006", + 'LT05_128016_20070515' = "2007", + 'LE07_128016_20070726' = "2007", + 'LE07_128016_20080626' = "2008", + 'LE07_128016_20080813' = "2008", + 'LE07_128016_20090715' = "2009", + 'LE07_128016_20090917' = "2009", + 'LE07_128016_20100819' = "2010", + 'LE07_128016_20110603' = "2010", + 'LT05_128016_20110713' = "2011", + 'LE07_128016_20120707' = "2012", + 'LE07_128016_20120723' = "2012", + 'LE07_128016_20130726' = "2013", + 'LE07_128016_20140611' = "2014", + 'LC08_128016_20140619' = "2014", + 'LC08_128016_20150606' = "2015", + 'LC08_128016_20150622' = "2015", + 'LE07_128016_20160803' = "2016", + 'LE07_128016_20160920' = "2016", + 'LC08_128016_20160608' = "2016", + 'LC08_128016_20160624' = "2016", + 'LC08_128016_20170611' = "2017", + 'LC08_128016_20170627' = "2017", + 'LE07_128016_20170705' = "2017", + 'LE07_128016_20170721' = "2017", + 'LC08_128016_20180614' = "2018", + 'LC08_128016_20180630' = "2018", + 'LE07_128016_20180521' = "2018", + 'LE07_128016_2018060' = "2018") + +year_ids <- c('layer.1'="1999", + 'layer.2'="2001", + 'layer.3'="2001", + 'layer.4'="2002", + 'layer.5'="2002", + 'layer.6'="2003", + 'layer.7'="2004", + 'layer.8'="2005", + 'layer.9'="2006", + 'layer.10'="2006", + 'layer.11'="2006", + 'layer.12'="2007", + 'layer.13'="2007", + 'layer.14'="2008", + 'layer.15'="2008", + 'layer.16'="2009", + 'layer.17'="2009", + 'layer.18'="2010", + 'layer.19'="2011", + 'layer.20'="2011", + 'layer.21'="2012", + 'layer.22'="2012", + 'layer.23'="2013", + 'layer.24'="2014", + 'layer.25'="2014", + 'layer.26'="2015", + 'layer.27'="2016", + 'layer.28'="2016", + 'layer.29'="2016", + 'layer.30'="2016", + 'layer.31'="2017", + 'layer.32'="2017", + 'layer.33'="2017", + 'layer.34'="2017", + 'layer.35'="2018", + 'layer.36'="2018", + 'layer.37'="2018", + 'layer.38'="2018") +``` + +```{r} +## Multipanel graph Summer 2000 - 2013 +lall_stack <- stack(LT05_127016_20000605, + LT05_127016_20010624, + LT05_127016_20020627, + LE07_127016_20050611, + LT05_127016_20060910, + LT05_127016_20070812, + LT05_127016_20090630, + LT05_127016_20100601, + LE07_127016_20120630, + LC08_127016_20130609) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +## Multipanel graph Augusts 2000 - 2018 +lall_stack <- stack(LT05_128016_19990728, + LT05_128016_20010701, + LT05_128016_20010717, + LT05_128016_20020704, + LT05_128016_20020720, + LE07_128016_20030715, + LT05_128016_20040725, + LE07_128016_20050704, + LE07_128016_20060520, + LE07_128016_20060723, + LT05_128016_20060731, + LT05_128016_20070515, + LE07_128016_20070726, + LE07_128016_20080626, + LE07_128016_20080813, + LE07_128016_20090715, + LE07_128016_20090917, + LE07_128016_20100819, + LE07_128016_20110603, + LT05_128016_20110713, + LE07_128016_20120707, + LE07_128016_20120723, + LE07_128016_20130726, + LE07_128016_20140611, + LC08_128016_20140619, + LC08_128016_20150606, + LC08_128016_20150622, + LE07_128016_20160803, + LE07_128016_20160920, + LC08_128016_20160608, + LC08_128016_20160624, + LC08_128016_20170611, + LC08_128016_20170627, + LE07_128016_20170705, + LE07_128016_20170721, + LC08_128016_20180614, + LC08_128016_20180630, + LE07_128016_20180521, + LE07_128016_20180606) + +lall_stack_df <- as.data.frame(lall_stack, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + facet_wrap(~ variable, labeller = as_labeller(year_id), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_timeseries_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2000_stack_ndvi_df <- as.data.frame(lall_stack[[1]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2000_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2000_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l1996_stack_ndvi_df <- as.data.frame(lall_stack[[9]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l1996_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_1996_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2016_stack_ndvi_df <- as.data.frame(lall_stack[[24]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_ndvi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(0,0.95), guide = guide_colorbar(title ="NDVI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/ndvi_2016_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + + +```{r} +# Calculate mean +lall_stack_mean <- calc(lall_stack, mean, na.rm=T) +# Calculate sd +lall_stack_sd <- calc(lall_stack, sd, na.rm=T) +lall_stack_emsi <- overlay(lall_stack, lall_stack_mean, lall_stack_sd, fun = function(r1, r2, r3) { return( (r1 - r2)/r3) }) +``` + +```{r} +# Rename variable layers back to year dates +year_ids <- c('layer.1'="2001", + 'layer.2'="2002", + 'layer.3'="2005", + 'layer.4'="2006", + 'layer.5'="2009", + 'layer.6'="2010", + 'layer.7'="2012", + 'layer.8'="2013") +``` + +```{r} +lall_stack_emsi_df <- as.data.frame(lall_stack_emsi, xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = lall_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + facet_wrap(~ variable, labeller = as_labeller(year_ids), ncol = 6) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_timeseries_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2010_stack_emsi_df <- as.data.frame(lall_stack_emsi[[8]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2010_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2010_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2007_stack_emsi_df <- as.data.frame(lall_stack_emsi[[6]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2007_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3.2,2.7), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2007_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + +```{r} +l2016_stack_emsi_df <- as.data.frame(lall_stack_emsi[[24]], xy = TRUE) %>% + melt(id.vars = c('x','y')) + +ggplot() + + geom_raster(data = l2016_stack_emsi_df, aes(x = x, y = y, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-3,3.2), guide = guide_colorbar(title ="EMSI")) + + theme(axis.title = element_blank(), + axis.text.x = element_blank(), + axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) + +ggsave('~/Documents/emsi_2016_yakutia.png', width = 12, height = 8, dpi = 300, bg = "transparent") +``` + + +```{r} +ggplot(lall_stack_emsi_df) + geom_histogram(aes(value, fill = value)) + + scale_fill_gradientn(colours=c("#6E462C", "#9C8448", "#CCCC66", "#9CAB68", "#306466"), limits=c(-2.5,2.7), guide = guide_colorbar(title ="EMSI")) + + ylab("Density") + xlab("EMSI") + ggtitle("June") + + facet_wrap(~variable, labeller = as_labeller(year_ids), ncol = 6) + + theme_bw() + + theme(axis.text.y = element_blank(), + axis.ticks = element_blank(), + rect = element_blank()) +``` \ No newline at end of file