Wednesday, 13 June 2018

javascript - Computing volumetric precipitation in Earth Engine



I am attempting to calculate volumetric precipitation (km^3) over a region in California, but the values seem too high to be plausible. The data suggests an average of 12 cm of precipitation for every point in the study area in February, 2014.


The code is below, and here is the link to the script


var area = ee.Geometry.Polygon([[[-118.654482361682, 37.1411680413078],
[-118.762206916723, 37.4566895272764],
[-119.259398430534, 37.7265422810242],
[-119.310905415871, 38.0449806440715],
[-119.629089966419, 38.1960150761288],

[-120.464922115121, 39.4371643262856],
[-120.105842345887, 39.5756453135706],
[-120.200754894698, 40.0567489586572],
[-121.254897840979, 40.5157526850281],
[-120.566012555589, 41.0835813539384],
[-120.200851104073, 41.0314786331859],
[-120.296983007049, 41.6665254634502],
[-120.187099478053, 41.9766005848439],
[-120.346617795513, 42.3774141613051],
[-120.836358420795, 42.3069364176645],

[-120.910857873804, 42.0594011597154],
[-120.563679730593, 41.9142972870239],
[-120.865943253041, 41.5290780147052],
[-121.024192398628, 41.590023082319],
[-121.118700296398, 41.4223695700792],
[-121.280739774272, 41.6293587801745],
[-121.589676536292, 41.6216515791448],
[-122.502956421333, 41.3093185577546],
[-122.448649964126, 41.1192597090914],
[-122.752298099071, 40.6899888128826],

[-122.693419681454, 40.574886307853],
[-123.068790630871, 40.3321184728131],
[-122.943612163357, 39.8016798173865],
[-122.668024701285, 39.5398636719596],
[-122.810965278147, 39.4087764148713],
[-122.711555463718, 39.256087663025],
[-123.024474210107, 39.319657188968],
[-123.096850519369, 39.0831079653768],
[-121.888502502495, 38.2885909322351],
[-121.936799168045, 37.8567862141554],

[-121.652246171612, 37.7636111215916],
[-121.407279611575, 37.1589020641969],
[-121.240316380584, 37.1587123798222],
[-121.234232158719, 36.9260955853916],
[-120.959661963312, 36.7198559190451],
[-121.007756475737, 36.5409587036978],
[-120.583862371187, 36.3328776238124],
[-120.66766576064, 36.1380952428648],
[-120.247896604, 35.8830368036773],
[-120.082411776132, 35.4949389115714],

[-119.350139889768, 34.8755948844077],
[-119.010002372588, 34.7768994366443],
[-118.24015921545, 35.2222860484529],
[-118.335924460093, 35.4188208481479],
[-117.985772699178, 35.6933520623051],
[-118.101414932332, 36.3460167748337],
[-118.436492979728, 37.0606299893495],
[-118.654482361682, 37.1411680413078]]])



var col = ee.ImageCollection("OREGONSTATE/PRISM/AN81m")
var t = col.filter(ee.Filter.calendarRange(2014, 2014, 'year')).filter(ee.Filter.calendarRange(2, 2, 'month')).select("ppt").filterBounds(area).sum()
var t2 = t.multiply(1e-3).multiply(ee.Image.pixelArea()).multiply(1e-9)
// This line does the following:
// (1) convert mm to m, (2) multiply by pixel area (given in m^2),(3) convert cubic m to cubic km

var scale = t2.projection().nominalScale()
var sumdict = t2.reduceRegion({
reducer:ee.Reducer.sum(),
geometry:area,

scale:scale})

print(sumdict) // = 19.6

// Average value of all precipitation over whole area = vol of preci / area
print(ee.Number(19.6).divide(area.area().multiply(1e-6)))
// = average of 12 cm precipitation for whole area. Seems high?

Map.addLayer(area)
Map.addLayer(t.clip(area),{"opacity":1,"bands":["ppt"],"min":15.191900234222413,"max":456.1631059837341,"gamma":1}, "prism P")


Answer



You should check this making a bar chart of multiple months/years with data existing. I think you are doing the job quite well. At least the trend of a lot of precipitation in the winter months and almost zero in June/July is very well visible.


To check your data for multiple months, you can use this double mapped function over your collection (first for years, then for months). Note that there exists a monthly dataset and a daily dataset. A added them both in this script and you can see that they are almost identical.


// Set years and month
var startYear = 2015;
var endYear = 2018;
var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(1,12);
// set the band(s) of interest
var bands = ['ppt'];


//////////// FOR THE DAILY DATASET \\\\\\\\\\\\\\\\\\\\\\
// loop over the years and months to get summed monthly images
var byMonth = years.map(function(y){
var yearCollection = Daily.filter(ee.Filter.calendarRange(y, y, 'year')).select(bands);
var byYear = ee.ImageCollection.fromImages(
months.map(function(m) {
var summedImage = yearCollection.filter(ee.Filter.calendarRange(m, m, 'month'))
.filterBounds(area).sum() //.set('month', m).set('year', y)
.set('system:time_start', ee.Date.fromYMD(y, m, 1, 'America/Los_Angeles'));

var t2 = summedImage.multiply(1e-3).multiply(ee.Image.pixelArea()).multiply(1e-9);
// This line does the following:
// (1) convert mm to m, (2) multiply by pixel area (given in m^2),(3) convert cubic m to cubic km
var sumDict = t2.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: area,
scale: t2.projection().nominalScale()});
// divide km^3 by area in km^2 and multiply by 1e5 to get precipitation in cm
var cmPrecPerMonth = ee.Number(sumDict.get('ppt')).divide(area.area().multiply(1e-6)).multiply(1e5);
return ee.Image(summedImage.setMulti(sumDict)).set('cmPrecPerMonth', cmPrecPerMonth);

}));
return byYear;
});

// rearrange the list of collections into a single image collection
var sizeList = byMonth.length();
var emptyCol = ee.ImageCollection(ee.Image());
function rearrange(current, previous){
previous = ee.ImageCollection(previous);
current = ee.ImageCollection(current);

return current.merge(previous);
}
// use iterate function and filter the empty one out
var outputDaily = ee.ImageCollection(byMonth.iterate(rearrange, emptyCol))
.filter(ee.Filter.listContains('system:band_names', 'ppt'))
.sort('system:time_start');
print('Monthly summed daily precipitation images', outputDaily);

print('Monthly summed daily precipitation images', ui.Chart.feature.byFeature(outputDaily, 'system:time_start', 'cmPrecPerMonth')
.setChartType('ColumnChart'));

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...