In [5]:
import shapely as sl
import shapely.geometry
import fiona
import ee
from utils.shapely_plot import draw
pylab.rcParams['figure.figsize'] = (20.0, 15.0)
ee.Initialize()
In [6]:
basins_au = [
None,
None,
ee.FeatureCollection('ft:1Dq_Q2JvvYkYO-kFX7L4E4Nzycwc50j9hfhSsBQJW'),
ee.FeatureCollection('ft:1778IyIZLZKSKVgko9X3aIV94E7jcm28uniyD6ycp'),
ee.FeatureCollection('ft:1WZ4Utbbatdl3vFVK7kTmAyHDyRjhMVfXeJeJTnBa'),
ee.FeatureCollection('ft:1rrk-yEOb8ILSolV_kSVD1qGxszHcy0cSL9UnUxIh'),
ee.FeatureCollection('ft:1-aMEhsi4usdxVUSSjKkJGC8pir3duCi_5oItnxtT'),
ee.FeatureCollection('ft:1YDeXF2LN8gDeJAOJTX0Kwp9QwV_-ZFI2llKilTGu'),
ee.FeatureCollection('ft:1YQ1qpXis4Z9z0NvKLdz-FjxFP5q2_fABi6aNSFn0')
];
In [7]:
# REGION1 Murray & Darling Catchment
# var region = { catchments: basins_au[3], main: 564, sub: basins_au[7], sub_min: 5640000, sub_max: 5649999, zoom: 6 }
# var region = { catchments: basins_au[3], main: 564, sub: basins_au[8], sub_min: 56400000, sub_max: 56499999, zoom: 6 }
# var region = { catchments: basins_au[4], main: 5641, sub: basins_au[8], sub_min: 56410000, sub_max: 56419999, zoom: 7 }
# var region = { catchments: basins_au[4], main: 5641, sub: basins_au[7], sub_min: 5641000, sub_max: 5641999, zoom: 7 }
region = { 'catchments': basins_au[6], 'main': 564105, 'sub': basins_au[8], 'sub_min': 56410500, 'sub_max': 56410599, 'zoom': 9 }
aoi_features = region['catchments'];
aoi = aoi_features.filter(ee.Filter.eq('PFAF_ID', region['main']));
not_aoi = aoi_features.filter(ee.Filter.neq('PFAF_ID', region['main']));
main_catchment_pfaf_min = region['sub_min'];
main_catchment_pfaf_max = region['sub_max'];
print('Area: ', aoi.geometry().area())
print('Pixels: ', aoi.geometry().area().getInfo() / (30.0*30.0))
center = aoi.geometry().centroid().getInfo()['coordinates']
print(center)
In [24]:
image = ee.Image('CGIAR/SRTM90_V4')
north = 37.0
south = 35.0
east = -111.0
west = -115.0
coords = [[west, north], [west, south], [east, south], [east, north], [west, north]]
url = image.visualize(min=[0], max=[3000]).getThumbUrl({'region': coords, 'format': 'png', 'size': '500'})
from IPython.display import Image
Image(url=url)
Out[24]:
In [18]:
image = aoi.reduceToImage(['PFAF_ID'], ee.Reducer.min())
?image.getThumbURL
url = image.getThumbURL({'dimensions':"800x600"})
print(url)
from IPython.display import Image
Image(filename=url)
In [5]:
Map.centerObject(aoi, region['zoom']);
var subcatchments = region.sub
.filter(ee.Filter.rangeContains('PFAF_ID', main_catchment_pfaf_min, main_catchment_pfaf_max))
//.limit(10);
var addBG = function() {
addToMapAsRaster(ee.FeatureCollection(ee.Geometry(Map.getBounds(true))), 'map (white)', 'ffffff', 0, 1, true, false);
addToMapAsRaster(aoi, 'aoi (black)', '000000', 0, 1, true, false);
}
var addSmallCatchments = function() {
addToMapAsRaster(not_aoi, 'not aoi', '000000,101010', 0, 0.5, true, true);
addToMapAsRaster(subcatchments, 'catchments', '101030,000000', 1, 0.9, false, false);
addToMapAsRaster(subcatchments, 'catchments (light)', 'a0a0ff,000000', 1, 0.9, false, false);
}
var LC8_BANDS = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'BQA'];
var LC7_BANDS = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'B8', 'B8', 'B7'];
var STD_NAMES = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pan', 'temp', 'BQA'];
var images = ee.ImageCollection('LANDSAT/LC8_L1T_TOA')
//.filterDate("2014-01-01", "2015-01-01")
.filterBounds(aoi)
//.filterBounds(Map.getBounds(true))
.select(LC8_BANDS, STD_NAMES);
var maxImageCount = 100;
// compute cloud cover
var addClouds = function() {
// total number of images
var count = images.select(0).count().clip(aoi);
Map.addLayer(count, {min:0, max:maxImageCount, palette:'cbc9e2,9e9ac8,6a51a3', opacity:0.9}, 'pixel count', false);
var chart = Chart.image.histogram(count, aoi, scale, 50);
chart = chart.setOptions({ title: 'Pixel Count' });
//chart.setChartType('PieChart');
print(chart)
// number of clean pixels using BQA
var bad = [61440,59424,57344,56320,53248,52256,52224,49184,49152,39936,31744,28590,26656]
var clearFn = function(img) { return img.select('BQA').eq(bad).reduce('max').not(); };
var clear = images.select('BQA').map(clearFn).sum().clip(aoi);
// Map.addLayer(clear, {min:0, max:maxImageCount, palette:'cbc9e2,9e9ac8,6a51a3', opacity:0.9}, 'clear', false)
var cloudRatio = ee.Image(1).subtract(clear.divide(count));
Map.addLayer(cloudRatio.mask(cloudRatio), {palette:'FFFFFF'}, 'clouds ratio', false)
//print(Chart.image.histogram(cloudRatio, aoi, 60, 80));
for(var th = 0.1; th <= 0.9; th += 0.1) {
Map.addLayer(cloudRatio.mask(cloudRatio.gte(th)).mask(cloudRatio.mask(cloudRatio.gte(th))), {min:0, max:1, palette:'FFFFFF'}, 'clouds > ' + Math.round(th*100) + '%', false);
}
// compute mean cloud ratio per sub-catchment
subcatchments = subcatchments.map(function(catchment){
var catchmentCloudRatio = cloudRatio.reduceRegion(ee.Reducer.mean(), catchment.geometry(), scale).get('constant');
return catchment.set('cloud_ratio', catchmentCloudRatio);
})
Map.addLayer(subcatchments.reduceToImage(['cloud_ratio'], ee.Reducer.first()), {palette: ['000000', 'ffffff']}, 'cloud ratio (mean)', false);
}
// compute NDWI min/max/mean, use cloud ratio to select percentiles
var waterIndices = function() {
var mean = images.reduce(ee.Reducer.intervalMean(20, 21));
//var ndwi = mean.normalizedDifference(['nir_mean', 'red_mean']) // NDVI
var ndwi = mean.normalizedDifference(['nir_mean', 'green_mean']) // NDWI
//var ndwi = mean.normalizedDifference(['swir1_mean', 'green_mean']) // MNDWI
var computeNDWI = function(catchment) {
var ndwiMin = ndwi.reduceRegion(ee.Reducer.min(), catchment.geometry(), 30).get('nd');
//var ndwiMax = ndwi.reduceRegion(ee.Reducer.max(), catchment.geometry(), 30).get('nd');
//var ndwiVar = ndwi.reduceRegion(ee.Reducer.variance(), catchment.geometry(), 30).get('nd');
return catchment
//.set('NDWI_var', ndwiVar)
.set('NDWI_min', ndwiMin)
//.set('NDWI_max', ndwiMax);
}
//Map.addLayer(ee.FeatureCollection([computeNDWI(ee.Feature(subcatchments.first()))]).reduceToImage(['NDWI_min'], ee.Reducer.first()), {min:-0.5, max:0.1, palette: ['000000', '0000ff']}, 'NDWI min', false);
// compute mean cloud ratio per sub-catchment
subcatchments = subcatchments.map(computeNDWI)
//Map.addLayer(subcatchments.reduceToImage(['NDWI_var'], ee.Reducer.first()), {min:0, max:0.1, palette: ['000000', 'ff0000']}, 'NDWI var', false);
Map.addLayer(subcatchments.reduceToImage(['NDWI_min'], ee.Reducer.first()), {min:-0.5, max:0.1, palette: ['000000', '0000ff']}, 'NDWI min', false);
//Map.addLayer(subcatchments.reduceToImage(['NDWI_max'], ee.Reducer.first()), {min:0.1, max:0.5, palette: ['000000', '00ff00']}, 'NDWI max', false);
//Map.addLayer(subcatchments.filter(ee.Filter.gt('NDWI_min', 0.2))
// .reduceToImage(['NDWI_min'], ee.Reducer.first()), {palette: ['ffffff']}, 'NDWI min > 0.2', false);
}
var addPercentiles = function() {
var images_sng = images.select('swir2', 'nir', 'green')
for(var i=10; i<100; i+=10) {
Map.addLayer(images_sng.reduce(ee.Reducer.intervalMean(i, i+1)).clip(aoi),
{min:0.05, max:0.5}, 'image(' + i + '%)', false)
}
}
//addBG()
//addPercentiles()
waterIndices()
//addClouds()
addSmallCatchments()
In [ ]:
sanFrancisco = ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8)
landsat8Toa = ee.ImageCollection('LANDSAT/LC8_L1T_32DAY_TOA').filterDate('2012-12-25', '2013-12-25').select('B[1-7]')
def getMean(img):
return img.reduceRegions(sanFrancisco, ee.Reducer.mean()
.forEachBand(img),200).makeArray(['B{0}'.format(x) for x in range(1,8)],'values')
ans = ee.FeatureCollection(landsat8Toa.map(getMean)).flatten().aggregate_array('.all').getInfo()
data = [x['properties']['values'] for x in ans]
x_ = range(0,7)
labels = ['B{0}'.format(x) for x in range(1,8)]
plt.plot(data)
plt.xticks(range(len(plt.xticks()[0])),[datetime.datetime.strptime(x['id'][:-2], "%Y%m%d").strftime('%b %Y') for x in ans])
plt.legend(plt.plot(data),['B{0}'.format(x) for x in range(1,8)])
plt.show()