-
Notifications
You must be signed in to change notification settings - Fork 0
/
ndvi.py
65 lines (50 loc) · 1.55 KB
/
ndvi.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
65
#!/usr/bin/python
from netCDF4 import Dataset
from mpl_toolkits.basemap import Basemap
from matplotlib import colors
import numpy as np
import matplotlib.pyplot as plt
#setting
clearthreshold = 296.0
#file
rawB03name = 'H08_B03_Indonesia_201706090500.nc'
rawB04name = 'H08_B04_Indonesia_201706090500.nc'
rawB13name = 'H08_B13_Indonesia_201706090500.nc'
#######################
#read raw data session#
#######################
dsetB03 = Dataset(rawB03name, mode = 'r')
dsetB04 = Dataset(rawB04name, mode = 'r')
dsetB13 = Dataset(rawB13name, mode = 'r')
#get lat and lon
lat = dsetB03.variables['latitude'][:]
lon = dsetB03.variables['longitude'][:]
#get data
B03 = dsetB03.variables['VS'][0]
B04 = dsetB04.variables['N1'][0]
B13 = dsetB13.variables['IR'][0]
###################
#calculate session#
###################
NDVI = (B04 - B03) / (B04 + B03)
NDVI[NDVI < 0.2] = np.nan
B13[B13 < clearthreshold] = np.nan
#NDVI = np.ma.masked_invalid(np.atleast_2d(NDVI))
NDVI = np.ma.masked_where(np.isnan(NDVI),NDVI)
#NDVI = np.ma.masked_where(np.isnan(B13),NDVI)
print NDVI
###############
#plotting time#
###############
m = Basemap(resolution='l', projection='merc', \
llcrnrlon=90, llcrnrlat=-20, urcrnrlon=149, urcrnrlat=20) #for mercator
#draw map coastline, etc
m.drawcoastlines(linewidth=0.75)
m.drawcountries(linewidth=0.75)
#make our lon and lat 2D
lon1, lat1 = np.meshgrid(lon, lat)
#for NDVI
plt.title('NDVI')
cs = m.pcolormesh(lon1, lat1, np.squeeze(NDVI), cmap='coolwarm', latlon=True)
cbar = m.colorbar(cs, location='bottom', pad="10%") #add color bar
plt.show()