<jupyter_start><jupyter_text>Computing NDVI over farmland Download your Landsat sceneUse the instructions at the beginning of `Lesson 1` to download another scene. This time, we're looking for some farm land. Take a guess at an area that is highly covered by farmable land, and download the __Red__ and __Infrared__ bands. Note: Record the date of the Landsat scene. We're going to use the date the Landsat scene was taken later, so write it down! Just like in the last lesson, we'll load up the raster data from our bands:<jupyter_code>import os

landsat_scene_name = "LC08_L1TP_042035_20171022_20171107_01_T1"

red_path = os.path.join("/home/hadoop/data/", 
                        "{}_B4.TIF".format(landsat_scene_name))
ir_path = os.path.join("/home/hadoop/data/", 
                       "{}_B5.TIF".format(landsat_scene_name))
import rasterio

with rasterio.open(red_path) as ds:
    red_data = ds.read()
    
with rasterio.open(ir_path) as ds:
    ir_data = ds.read()<jupyter_output><empty_output><jupyter_text>Now let's use some of `numpy`'s computational capabilites to do some band math. We'll be computing the [Normalized difference vegetation index](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index), which computes the ratio of the difference and the sum of the Red and Near Infrared bands of Landsat to create a value that has proven to be an efficient indicator of vegetation health. The formula for NDVI is:![NDVI](img/ndvi.png)The range for this function is `-1.0` to `1.0`. The greater the number, the greater the indication of healthy vegetation.To encode this in `numpy` operations, first we will change our values to floating point so that the result of the operation is also a floating point value:<jupyter_code>import numpy as np

r = red_data[0].astype(float)
ir = ir_data[0].astype(float)<jupyter_output><empty_output><jupyter_text>Then we utilize numpy's overriding of the division, subtraction and addition operators to write code that looks a whole lot like the original formula:<jupyter_code>ndvi = (ir - r) / (ir + r)<jupyter_output><empty_output><jupyter_text>You may notice a warning that occurs, which happens when the division hits a `0` value for the `ir` and `r` data. This happens because `0` is used as a `nodata` value for Landsat. To get around this, we can use `numpy` syntax to set every 0 to `nan` (Not A Number):<jupyter_code>r = red_data[0].astype(float)
r[r == 0.0] = np.nan
ir = ir_data[0].astype(float)
ir[ir == 0.0] = np.nan<jupyter_output><empty_output><jupyter_text>Now, the division does not throw this warning:<jupyter_code>ndvi = (ir - r) / (ir + r)<jupyter_output><empty_output><jupyter_text>We'll be computing NDVI a number of times, so let's refactor the above code into a function:<jupyter_code>def compute_ndvi(r, ir):
    rf = r.astype(float)
    rf[r == 0] = np.nan
    irf = ir.astype(float)
    irf[ir == 0] = np.nan
    return (irf - rf) / (irf + rf)
ndvi = compute_ndvi(red_data[0], ir_data[0])<jupyter_output><empty_output><jupyter_text>Statistics of the NDVIHere we inspect statistics of the NDVI values:<jupyter_code>flattened = np.ravel(ndvi[~np.isnan(ndvi)])

print("MEAN: {}".format(np.mean(flattened)))
print("MEDIAN: {}".format(np.median(flattened)))
print("MIN: {}".format(np.min(flattened)))
print("MAX: {}".format(np.max(flattened)))
import matplotlib.pyplot as plt
plt.hist(flattened, bins='auto')  # arguments are passed to np.histogram
plt.title("Histogram with 'auto' bins")
plt.show()<jupyter_output><empty_output><jupyter_text>The rendering of the histogram should give you an idea of the normal range of NDVI values.  Clipping out a subset of fieldsWe're now going to use [geojson.io](http://geojson.io) to find a subset of the data to work with, that contains fields for us to analyze.Let's get the bounds of the image, as well as the CRS from the dataset:<jupyter_code>with rasterio.open(red_path) as ds:
    bounds = ds.bounds
    crs = ds.crs<jupyter_output><empty_output><jupyter_text>Now, we can create a `shapely` polygon out of the bounds, and use the `mapping` method to turn it into a GeoJSON `dict`:<jupyter_code>from shapely.geometry import Polygon

bounds_poly = Polygon([(bounds.left, bounds.bottom),
                       (bounds.left, bounds.top),
                       (bounds.right, bounds.top),
                       (bounds.right, bounds.bottom),
                       (bounds.left, bounds.bottom)])

bounds_poly
from shapely.geometry import mapping

m = mapping(bounds_poly)
m<jupyter_output><empty_output><jupyter_text>You'll notice the coordinates don't look anything like the lat/long coordinates from the county in Lesson 1. Again, we have to mind our projections. These coordinates are in the image's CRS, so let's reproject them to lat/long:<jupyter_code># Reproject to LatLng
import pyproj
from musa import reproject_geom # From Lesson 1

image_proj = pyproj.Proj(crs)
lat_lng_proj = pyproj.Proj(init="epsg:4326")
ll_bounds_poly = reproject_geom(bounds_poly, image_proj, lat_lng_proj)
ll_bounds_poly<jupyter_output><empty_output><jupyter_text>Notice I imported the `reproject_geom` method from a `musa`. This is a module that lives in the workshop repository, in which I've put some methods that either we have already gone over or aren't valuable to the lesson at hand.Now, we can print out the GeoJSON and load it up in geojson.io:<jupyter_code># Dump to GeoJSON
import json

geojson = json.dumps(mapping(ll_bounds_poly), indent=4)
print(geojson)<jupyter_output><empty_output><jupyter_text>Copy the GeoJSON and paste it into [geojson.io](http://geojson.io). It should look something like this:![geojson.io](img/geojson-io-1.png)Zoom into a part of the land covered by the polygon that seems likely to be covered by farmland. It's useful to hit the "Satellite" basemap toggle on the bottom left. Once you've zoomed to your area, use the trash icon to delete the bounds polygon:![geojson.io](img/geojson-io-2.png)Now, use the polygon icon to draw a polygon around some fields. Try to make the area big enough that it covers a lot of farmland, but not so big that we might as well not be subsetting: ![geojson.io](img/geojson-io-3.png) Copy out all the text in the GeoJSON editor and paste it to replace `` below:<jupyter_code># Load in GeoJSON

crop_area_geojson = """
<YOUR TEXT HERE>
"""<jupyter_output><empty_output><jupyter_text>Now, we can grab the geometry from the FeatureCollection of the GeoJSON, reproject to the Landsat scene's CRS, and use it as a mask to  crop our NDVI:<jupyter_code>from shapely.geometry import shape

crop_area_json = json.loads(crop_area_geojson)['features'][0]['geometry']
crop_area_ll = shape(crop_area_json)
crop_area_ll
crop_area = reproject_geom(crop_area_ll, lat_lng_proj, image_proj)
from rasterio.mask import mask

with rasterio.open(red_path) as ds:
    (cropped_red_data, cropped_transform) = mask(ds, [mapping(crop_area)], crop=True)
    crs = ds.crs
    out_meta = ds.meta.copy()

# Save a GeoTiff
with rasterio.open(ir_path) as ds:
    (cropped_ir_data, _) = mask(ds, [mapping(crop_area)], crop=True)
cropped_ndvi = compute_ndvi(cropped_red_data[0], cropped_ir_data[0])<jupyter_output><empty_output><jupyter_text>Because we targeted farm land, one might think the stats for NDVI will shift towards `1.0`. However, because many fields might either not be in growing season or have already been harvested, you might find the numbers don't change that much or may even skew downward.<jupyter_code>cropped_flattened = np.ravel(cropped_ndvi[~np.isnan(cropped_ndvi)])
print("MEAN: {}".format(np.mean(cropped_flattened)))
print("MEDIAN: {}".format(np.median(cropped_flattened)))
print("MIN: {}".format(np.min(cropped_flattened)))
print("MAX: {}".format(np.max(cropped_flattened)))
plt.hist(cropped_flattened, bins='auto')  # arguments are passed to np.histogram
plt.title("Histogram with 'auto' bins")
plt.show()<jupyter_output><empty_output><jupyter_text>Viewing NDVILet's take a look at the NDVI. NDVI is often viewed with a color palette that shows more intense green where there is healthy vegitation. Here we define a matplotlib color map based on blues and greens:<jupyter_code>from matplotlib.colors import LinearSegmentedColormap

# Divergent Blue to Beige to Green colormap
cmap = LinearSegmentedColormap.from_list(
  'ndvi', ['blue', 'beige', 'green'], 20)<jupyter_output><empty_output><jupyter_text>We can then show the image rendered with this color map:<jupyter_code>from musa import show_image # From Lesson 1

show_image(cropped_ndvi, cmap=cmap)<jupyter_output><empty_output><jupyter_text>Already you might be able to see distinct fields that are either harvested (having low NDVI values) or currently growing (having high NDVI values).  Writing a GeoTiffWith rasterio, you can write out data to a GeoTIFF. Here we write out our NDVI data, so that we can for example view it in [QGIS](https://qgis.org/en/site/) or process it further with [GDAL](http://www.gdal.org/):<jupyter_code>out_image = np.expand_dims(ndvi, axis=0)

out_meta.update({"driver": "GTiff",
                 "dtype": "float64",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": cropped_transform})

with rasterio.open("/home/hadoop/data/croptest.tif", "w", **out_meta) as dest:
    dest.write(out_image)<jupyter_output><empty_output><jupyter_text>Viewing with GeoPySpark and GeoNotebookAnother option to view is to use [GeoPySpark](https://geopyspark.readthedocs.io/en/latest/) put the data on this map __→__We haven't utilized this map and GeoPySpark and [GeoNotebook](https://github.com/OpenGeoscience/geonotebook) for a couple of reasons, mainly because these topics are out of scope for the workshop. However, let's use some functionality I put into the `musa` code to render our NDVI onto the map.First we take the bounds of image, center the map on it, and call the `map_ndvi` method to render it to the map. This might take a few minutes. For a peak under the hood, check out the python code in the `musa` directory of this repo.<jupyter_code>from rasterio.transform import array_bounds
from rasterio.coords import BoundingBox

b = array_bounds(cropped_red_data.shape[1], cropped_red_data.shape[2], cropped_transform)
bounds = BoundingBox(b[0], b[1], b[2], b[3])
from musa import map_ndvi
center = crop_area_ll.centroid
M.set_center(center.x, center.y, z = 10)<jupyter_output><empty_output><jupyter_text>If your data isn't small or your machine is underpowered, this next step can take a bit. If it's taking a while, or you are interested, go to http://locahost:4040 and check on the spark job. Each tile is served out of an Apache Spark RDD from an in-memory tile server.<jupyter_code># Remove any previously placed on layers
for l in M.layers:
    M.remove_layer(l)

map_ndvi(M, cropped_ndvi, bounds, crs)<jupyter_output><empty_output><jupyter_text>Finding the NDVI values of a single fieldWe'll use the GeoNotebook functionality for drawing a polygon on the  map, and use it to single out an individual field in our data.Click the draw polygon button in the toolbar of this notebook. It looks like this:![Draw Polygon](img/geonotebook-draw.png)Then draw a polygon around a field that looks particularly healthy.Now we can grab the polygon using the GeoNotebook annotation layer:<jupyter_code>aoi = M.layers.annotation.polygons[0]
aoi<jupyter_output><empty_output><jupyter_text>Rasterio has a `features.geom_mask` function that will create a mask of raster data based on a geometry you pass in. We can use this to mask out the part of the image that is not covered by our polygon, and then show a histogram of values that represent our field:<jupyter_code>from rasterio import features

r_aoi = reproject_geom(aoi, lat_lng_proj, image_proj)
geom_mask = features.geometry_mask(
        [r_aoi],
        out_shape=cropped_ndvi.shape,
        transform=cropped_transform
    )
masked = np.ma.array(data=cropped_ndvi, mask=geom_mask)
np.mean(masked)
plt.hist(np.ravel(cropped_ndvi[~np.isnan(cropped_ndvi)]), bins='auto')  # arguments are passed to np.histogram
plt.title("Histogram with 'auto' bins")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09493307054282664, 'fraction_numerical': 0.0070624948673729164, 'mean_word_length': 5.089, 'pattern_counts': {'":': 5, '<': 86, '<?xml version=': 0, '>': 86, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '470480', 'license': 'Apache-2.0', 'n_tokens_mistral': 4115, 'n_tokens_neox': 3918, 'n_words': 1361, 'path': 'notebooks/Lesson 2 - Computing NDVI over an area.ipynb', 'repo_name': 'zhaoanbei/raster-processing-with-python-musa-2018', 'size': '12177'}