-
Notifications
You must be signed in to change notification settings - Fork 49
/
idw_interpolation.py
64 lines (53 loc) · 1.89 KB
/
idw_interpolation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import ee
import geemap
# Create a map centered at (lat, lon).
Map = geemap.Map(center=[40, -100], zoom=4)
def sampling(sample):
lat = sample.get('latitude')
lon = sample.get('longitude')
ch4 = sample.get('ch4')
return ee.Feature(ee.Geometry.Point([lon, lat]), {'ch4': ch4})
# Import two weeks of S5P methane and composite by mean.
ch4 = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_CH4') \
.select('CH4_column_volume_mixing_ratio_dry_air') \
.filterDate('2019-08-01', '2019-08-15') \
.mean() \
.rename('ch4')
# Define an area to perform interpolation over.
aoi = ee.Geometry.Polygon(
[[[-95.68487605978851, 43.09844605027055],
[-95.68487605978851, 37.39358590079781],
[-87.96148738791351, 37.39358590079781],
[-87.96148738791351, 43.09844605027055]]], {}, False)
# Sample the methane composite to generate a FeatureCollection.
samples = ch4.addBands(ee.Image.pixelLonLat()) \
.sample(**{'region': aoi, 'numPixels': 1500,
'scale':1000, 'projection': 'EPSG:4326'}) \
.map(sampling)
# Combine mean and standard deviation reducers for efficiency.
combinedReducer = ee.Reducer.mean().combine(**{
'reducer2': ee.Reducer.stdDev(),
'sharedInputs': True})
# Estimate global mean and standard deviation from the points.
stats = samples.reduceColumns(**{
'reducer': combinedReducer,
'selectors': ['ch4']})
# Do the interpolation, valid to 70 kilometers.
interpolated = samples.inverseDistance(**{
'range': 7e4,
'propertyName': 'ch4',
'mean': stats.get('mean'),
'stdDev': stats.get('stdDev'),
'gamma': 0.3})
# Define visualization arguments.
band_viz = {
'min': 1800,
'max': 1900,
'palette': ['0D0887', '5B02A3', '9A179B', 'CB4678',
'EB7852', 'FBB32F', 'F0F921']}
# Display to map.
# Map.centerObject(ee.FeatureCollection(aoi), 7)
Map.addLayer(ch4, band_viz, 'CH4')
# Map.addLayer(interpolated, band_viz, 'CH4 Interpolated')
# Display the map.
Map