# Set matplotlib plots inline
%matplotlib inline

# Import libraries
import urllib
import shutil
import os
import glob
import numpy as np
from netCDF4 import Dataset

# Set input directory, and change working directory
in_dir = input("Main input directory: ")                     # IMPORTANT: Need to enter your working directory
os.chdir(in_dir)

# Create and set output directory
out_dir = os.path.normpath((os.path.split(in_dir)[0] + os.sep + 'output/' ))+ '\\'
if not os.path.exists(out_dir):
    os.makedirs(out_dir)

# Search for, open, and read text file containing links to AppEEARS output files
download_file = open(glob.glob('**.txt')[0],'r')
download_link = download_file.readlines()

# Loop through text file and download all AppEEARS outputs
for i in range(len(download_link)):  
    file_name = download_link[i].split('/')[-1].strip()

    # Download the file
    response = urllib.request.urlopen(download_link[i])
    out_file = open(file_name, 'wb')
    shutil.copyfileobj(response, out_file)
    print('Downloaded file: '+ file_name + ' (' + str(i+1) + " of " + str(len(download_link)) + ')')
    out_file.close()

# Once you have finished downloading the files, remove variables not needed for the remainder of the tutorial
del download_link, file_name, download_file

Downloaded file: MCD12Q1.006_500m_aid0001.nc (1 of 2)
Downloaded file: MOD13A2.006_1km_aid0001.nc (2 of 2)

# Now that we have files downloaded, search for downloaded files
file_list = glob.glob('**.nc')
file_list

['MCD12Q1.051_aid0001.nc', 'MOD13A2.006_aid0001.nc']

# Read file in, starting with MOD13A2 version 6
file_in = Dataset(file_list[1],'r', format = 'NETCDF4')

# Print a list of variables in file
list(file_in.variables)

['crs', 'time', 'lat', 'lon', '_1_km_16_days_EVI', '_1_km_16_days_VI_Quality']

# Grab variable and variable metadata
v6_info = file_in.variables['_1_km_16_days_EVI']
print(v6_info)

# File dimensions
file_in.dimensions.values()

<class 'netCDF4._netCDF4.Variable'>
int16 _1_km_16_days_EVI(time, lat, lon)
    _FillValue: -3000
    coordinates: time lat lon
    grid_mapping: crs
    valid_min: -2000
    valid_max: 10000
    long_name: 1 km 16 days EVI
    units: EVI
    scale_factor_err: 0.0
    add_offset_err: 0.0
    calibrated_nt: 5
    scale_factor: 0.0001
    add_offset: 0.0
unlimited dimensions: 
current shape = (7, 3600, 4200)
filling on
dict_values([<class 'netCDF4._netCDF4.Dimension'>: name = 'time', size = 7, <class 
'netCDF4._netCDF4.Dimension'>: name = 'lat', size = 3600, <class 'netCDF4._netCDF4.Dimension'>: name = 
'lon', size = 4200])

from netCDF4 import num2date

times = file_in.variables["time"]                            # import the time variable 
dates = num2date(times[:], times.units)                      # get the time info
dates = [date.strftime('%Y-%m-%d') for date in dates]        # get the list of datetime
print(dates)

['2005-07-12', '2005-07-28', '2005-08-13', '2005-08-29', '2005-09-14', '2005-09-30']

# import the EVI and VI Quality datasets
v6_EVI = file_in.variables['_1_km_16_days_EVI'][1:]          # select the EVI layers from July to September 
print(v6_EVI.shape)
v6_QA = file_in.variables['_1_km_16_days_VI_Quality'][1:]   # select the EVI quality layers from July to September
print(v6_QA.shape)

(6, 3600, 4200)
(6, 3600, 4200)

# Set lat and lon min/max for EVI data
lat_EVI_min = file_in.variables['lat'][:].min()
lat_EVI_max = file_in.variables['lat'][:].max()
lon_EVI_min = file_in.variables['lon'][:].min()
lon_EVI_max = file_in.variables['lon'][:].max()

# Search for look up table 
lut = glob.glob('**lookup.csv')
lut

['MOD13A2-006-1-km-16-days-VI-Quality-lookup.csv']

import pandas as pd
# Read in the lut
v6_QA_lut = pd.read_csv(lut[0])
v6_QA_lut

# Include good quality based on MODLAND
modland = ['VI produced, good quality', 'VI produced, but check other QA']
v6_QA_lut = v6_QA_lut[v6_QA_lut['MODLAND'].isin(modland)]

# Exclude lower quality VI usefulness
VIU =["Lowest quality","Quality so low that it is not useful","L1B data faulty","Not useful for any other reason\/not processed"]
v6_QA_lut = v6_QA_lut[~v6_QA_lut['VI Usefulness'].isin(VIU)]

# Include low or average aerosol
AQ = ['Low','Average']
v6_QA_lut = v6_QA_lut[v6_QA_lut['Aerosol Quantity'].isin(AQ)]

# Include where adjacent cloud, mixed clouds, or possible shadow were not detected
v6_QA_lut = v6_QA_lut[v6_QA_lut['Adjacent cloud detected'] == 'No' ]
v6_QA_lut = v6_QA_lut[v6_QA_lut['Mixed Clouds'] == 'No' ]
v6_QA_lut = v6_QA_lut[v6_QA_lut['Possible shadow'] == 'No' ]
v6_QA_lut

# Print list of possible QA values based on parameters above
v6_QA_mask = list(v6_QA_lut['Value'])
del v6_QA_lut

v6_QA_mask

[2112,
 2116,
 2177,
 2181,
 2185,
 4160,
 4164,
 4225,
 4229,
 4233,
 6208,
 6212,
 6273,
 6277,
 6281]

# Apply QA mask to the EVI data 
v6_EVI_masked = np.ma.MaskedArray(v6_EVI[0], np.in1d(v6_QA[0], v6_QA_mask, invert = True))

import matplotlib.pyplot as plt

# Visualize a basic plot
plt.figure(figsize = (10,10))                # Set the figure size
plt.axis('off')                              # Remove the axes' values

# Plot the original data, using a colormap and setting a custom linear stretch based on the min/max EVI values
plt.imshow(v6_EVI[0], vmin = -0.2 , vmax = 1, cmap = 'YlGn');        # Adding ';' to the end of the call suppresses the text output from showing in the notebook 

# Set the figure size
fig = plt.figure(figsize = (10,10))

# Remove the axes' values
plt.axis('off')

# Make a figure title
fig.suptitle('MODIS Version 6 Enhanced Vegetation Index (EVI) Quality Filtered Data', fontsize = 12, fontweight = 'bold')

# Make a figure subtitle
ax = fig.add_subplot(111)
fig.subplots_adjust(top = 2.2)
ax.set_title('Amazon Rainforest: 07-12-2005')
ax1 = plt.gca()

# Plot the masked data
im = ax1.imshow(v6_EVI_masked, vmin = -0.2, vmax = 1, cmap = 'YlGn');

# Add a colormap legend
plt.colorbar(im, orientation='horizontal', fraction = 0.047, pad = 0.02, label = 'EVI').outline.set_visible(True)

# Set up file name and export to png file
import warnings
warnings.filterwarnings("ignore") 
fig.savefig('{}{}_evi.png'.format(out_dir, (file_list[1].split('.nc'))[0]))

# Use MCD12Q1.006 to define land cover types
file_list

['MCD12Q1.006_500m_aid0001.nc', 'MOD13A2.006_1km_aid0001.nc']

# set to 'MCD12Q1.006_aid0001.nc' and open
lct_file = Dataset(file_list[0],'r', format = 'NETCDF4')

# print a list of variables in file
list(lct_file.variables)

['crs', 'time', 'lat', 'lon', 'LC_Type1', 'QC']

# Show the variable metadata
lct_file.variables['LC_Type1']

<class 'netCDF4._netCDF4.Variable'>
int16 LC_Type1(time, lat, lon)
    _FillValue: 255
    coordinates: time lat lon
    grid_mapping: crs
    valid_min: 1
    valid_max: 17
    long_name: Land_Cover_Type_1
    Evergreen_Needleleaf_Forests: 1
    Evergreen_Broadleaf_Forests: 2
    Deciduous_Needleleaf_Forests: 3
    Deciduous_Broadleaf_Forests: 4
    Mixed_Forests: 5
    Closed_Shrublands: 6
    Open_Shrublands: 7
    Woody_Savannas: 8
    Savannas: 9
    Grasslands: 10
    Permanent_Wetlands: 11
    Croplands: 12
    Urban_and_Built_up_Lands: 13
    Cropland_Natural_Vegetation_Mosaics: 14
    Permanent_Snow_and_Ice: 15
    Barren: 16
    Water_Bodies: 17
    Unclassified: -1
unlimited dimensions: 
current shape = (1, 7200, 8400)
filling on

# Open the Land_Cover_Type_1 dataset
lct = lct_file.variables['LC_Type1'][:][0]

# Make a list of the forest land cover types outlined above
lct_forest = [1,2,3,4,5]

# Mask the dataset to only include forest LCTs
lct_masked_forest = np.ma.MaskedArray(lct, np.in1d(lct, lct_forest, invert = True))

import scipy.ndimage

# Resample by a factor of 2 with nearest neighbor interpolation
resampled_lct = scipy.ndimage.zoom(lct, 0.5, order=0)

# Next, exclude non-Forest land covers from the EVI data
v6_EVI_masked_forest = np.ma.MaskedArray(v6_EVI_masked, np.in1d(resampled_lct, lct_forest, invert = True))

del lct

# Show before/after quality screening
figure, axes = plt.subplots(nrows=1, ncols=3, figsize=(40, 10), sharey=True)
figure.suptitle('MODIS Version 6 Enhanced Vegetation Index (EVI): 07-12-2005', fontsize=32, fontweight='bold')

# Turn off axes' values
axes[0].axis('off'), axes[1].axis('off'), axes[2].axis('off')

# Plot the original data
axes[0].set_title('Original Data', fontsize=20, fontweight='bold')
axes[0].imshow(v6_EVI[0], vmin = -0.2, vmax = 1.0, cmap = 'YlGn')

# Plot the quality filtered data
axes[1].set_title('Quality Filtered Data', fontsize=20, fontweight='bold')
axes[1].imshow(v6_EVI_masked, vmin = -0.2, vmax = 1.0, cmap = 'YlGn')

# Plot the forest-only data
axes[2].set_title('Quality Filtered Forest Data', fontsize=20, fontweight='bold')
axes[2].imshow(v6_EVI_masked_forest, vmin = -0.2, vmax = 1.0, cmap = 'YlGn');

from matplotlib.colors import LinearSegmentedColormap

# Define colormap
lct_cmap = ((0.149,0.451,0),(0.220,0.659,0),(0.298,0.902,0),(0.639,1,0.451),(0.266,0.537,0.439))
cmap1 = LinearSegmentedColormap.from_list("LCT", lct_cmap, N=5, gamma=1.0)

del lct_cmap

# Remove Matplotlib deprecation warnings when importing Basemap
import warnings
warnings.simplefilter("ignore")

    projection = the desired projection (EPSG: 4326 here)
    ellps = the desired datum (WGS84 here)
    llcrnrlat= lower left bounding box latitude
    urcrnrlat= upper right bounding box latitude
    llcrnrlon= lower left bounding box longitude
    urcrnrlon= upper right bounding box longitude
    resolution= resolution of ancillary boundary layers in basemap plot, 'i' here = intermediate resolution

# Import Basemap library
from mpl_toolkits.basemap import Basemap ;

# Set up the figure
fig = plt.figure(figsize=(8,8))

# Set the extent values
x0, y0 = (lon_EVI_min, lat_EVI_min)
x1, y1 = (lon_EVI_max, lat_EVI_max)

# Create Basemap
m = Basemap(projection='cyl', ellps = 'WGS84',llcrnrlat=-20,urcrnrlat=13,llcrnrlon=-83,urcrnrlon=-34,resolution='l');

# Draw coastlines and country boundaries, zorder is used to define the order of elements
m.drawcoastlines(zorder = 10), m.drawcountries(zorder = 10)
m.drawmapboundary(fill_color='lightblue'), m.fillcontinents(color='white')

# Plot the LCT data
plt.imshow(lct_masked_forest, vmin = 1 , vmax = 5, cmap = cmap1, zorder =5, origin = 'upper', extent= (x0, x1, y0, y1))

# Add a title
plt.title("Forest Land Cover Types of the Amazon Basin (2005)")
plt.show()

warnings.simplefilter("ignore")
lct_file.close()
del lct_masked_forest

# Call np.mean() and print results
print('Version 6 all quality mean EVI: ' + str(float(format((np.mean(v6_EVI[0])), '.4f'))))
print('Version 6 quality filtered mean EVI: ' + str(float(format((np.mean(v6_EVI_masked)), '.4f'))))
print('Version 6 quality filtered mean (forests) EVI: ' + str(float(format((np.mean(v6_EVI_masked_forest)), '.4f'))))

Version 6 all quality mean EVI: 0.4078
Version 6 quality filtered mean EVI: 0.4072
Version 6 quality filtered mean (forests) EVI: 0.5003

# Make a copy of the original NetCDF-4 file
v6_EVI_final = v6_EVI

# Loop through all timesteps (observations) in the file, mask out poor quality and exclude non-forest pixels
for i in range(len(dates)):

    # Next, apply QA mask to the EVI data
    v6_EVI_masked = np.ma.MaskedArray(v6_EVI[i], np.in1d(v6_QA[i], v6_QA_mask, invert = True));

    # Next, # Next, exclude npn-Forest land covers from the EVI data
    v6_EVI_final[i] = np.ma.MaskedArray(v6_EVI_masked, np.in1d(resampled_lct, lct_forest, invert = True));
del v6_EVI, v6_EVI_masked, v6_QA, resampled_lct

# Set up subplots with row and column numbers
figure, axes = plt.subplots(nrows=int(len(dates)/3), ncols=3, figsize=(25, 10), sharey=True, sharex=True)
# Add a title for the series of plots
figure.suptitle('MODIS Version 6 Enhanced Vegetation Index (EVI): 2005-07-12 to 2005-09-30', fontsize=32, fontweight='bold')

for i in range(int(len(dates)/3+1)):
    for j in range(0,3):
        try:
            axes[i][j].set_title(dates[3*i+j],fontsize=20,fontweight='bold')
            # Plot each observation
            axes[i][j].imshow(v6_EVI_final[3*i+j],vmin = -0.2, vmax = 1.0, cmap = 'YlGn',zorder = 5, origin = 'upper', 
                              extent= (x0, x1, y0, y1))
            # Set up the basemap for each subplot
            m2 =Basemap(projection='cyl',llcrnrlat=-20,urcrnrlat=13,llcrnrlon=-83,urcrnrlon=-34,resolution='l',
                        ellps ='WGS84',ax=axes[i][j])
            m2.drawcoastlines(zorder=10)
            m2.drawcountries(zorder=10)
            m2.drawmapboundary(fill_color='lightblue'),m2.fillcontinents(color='white');
        except IndexError:
            pass

# Create output file
v6_file_out = Dataset(out_dir+file_list[1].split('.nc')[0] + '_masked.nc','w', format = 'NETCDF4')

#Copy dimensions of the input file
print(file_in.dimensions.items())
v6_file_out.createDimension('time', len(dates))
v6_file_out.createDimension('lat', len(file_in.dimensions['lat']))
v6_file_out.createDimension('lon', len(file_in.dimensions['lon']))

# Copy variables from the input file
v6_variables = list(file_in.variables)

# loop through each variable
for j in range(len(v6_variables)):
    # Copy the input variable 
    in_var = file_in.variables[v6_variables[j]]
    # Create the output variable
    out_var = v6_file_out.createVariable(v6_variables[j], in_var.datatype, in_var.dimensions)
    # Copy variable attributes
    out_var.setncatts({k: in_var.getncattr(k) for k in in_var.ncattrs()})
    # replace the original EVI data with our quality and forest masked data
    if v6_variables[j] == '_1_km_16_days_EVI':
        out_var[:] = v6_EVI_final[:]
    elif v6_variables[j] == 'time':
            out_var[:] = in_var[1:]
    elif v6_variables[j] =='_1_km_16_days_VI_Quality':
        out_var[:] = in_var[1:]
    else: 
        out_var[:] = in_var[:]
# close the input and output files
v6_file_out.close()
file_in.close()