In this tutorial, we will walk through how to remove parts of a raster based on pixel values using a mask we create. A mask raster layer contains pixels that are not used in the analysis. In Python, these pixels are assigned a value of nan (not a number).
In [1]:
import numpy as np
import gdal
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
In [2]:
# Import TEAK aspect raster generated in previous lesson
aspectNS_filename = 'teak_nsAspect.tif'
aspectNS_dataset = gdal.Open(aspectNS_filename)
#Display the dataset dimensions:
rows_aspNS = aspectNS_dataset.RasterYSize; print('x:',rows_aspNS)
cols_aspNS = aspectNS_dataset.RasterXSize; print('y:',cols_aspNS)
print('bands:',aspectNS_dataset.RasterCount)
print('driver:',aspectNS_dataset.GetDriver().LongName)
print('geotransform:',aspectNS_dataset.GetGeoTransform())
print('projection:',aspectNS_dataset.GetProjection())
aspectNS_mapinfo = aspectNS_dataset.GetGeoTransform()
xMin = aspectNS_mapinfo[0]
yMax = aspectNS_mapinfo[3]
xMax = xMin + aspectNS_dataset.RasterXSize/aspectNS_mapinfo[1] #divide by pixel width
yMin = yMax - aspectNS_dataset.RasterYSize/aspectNS_mapinfo[5] #divide by pixel height (note sign +/-)
aspectNS_ext = (xMin,xMax,yMin,yMax)
print('TEAK Aspect NS Raster Extent:',aspectNS_ext)
#Convert CHM raster to NumPy array
aspectNS_raster = aspectNS_dataset.GetRasterBand(1)
noDataVal = aspectNS_raster.GetNoDataValue()
print('no data value:',noDataVal)
scaleFactor = aspectNS_raster.GetScale()
print('scale factor:',scaleFactor)
aspectNS_stats = aspectNS_raster.GetStatistics(True,True)
print('CHM Statistics: Minimum=%.2f, Maximum=%.2f, Mean=%.3f, StDev=%.3f' %
(aspectNS_stats[0], aspectNS_stats[1], aspectNS_stats[2], aspectNS_stats[3]))
aspectNS_array = aspectNS_dataset.GetRasterBand(1).ReadAsArray(0,0,cols_aspNS,rows_aspNS).astype(np.float)
print(aspectNS_array)
In [3]:
# Plot classified aspect (N-S) array
from matplotlib import colors
aspectNS_array = aspectNS_array[::-1] #?? do raster arrays need to be inversed to look right as geotif?
fig, ax = plt.subplots()
cmapNS = colors.ListedColormap(['white','blue','red'])
plt.imshow(aspectNS_array,extent=aspectNS_ext,cmap=cmapNS)
plt.title('TEAK North & South Facing Slopes')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
# Create custom legend to label N & S
import matplotlib.patches as mpatches
white_box = mpatches.Patch(color='white',edgecolor='red',label='East, West, or NaN')
blue_box = mpatches.Patch(color='blue', label='North')
red_box = mpatches.Patch(color='red', label='South')
ax.legend(handles=[white_box,blue_box,red_box],handlelength=0.7,bbox_to_anchor=(1.05, 0.45),
loc='lower left', borderaxespad=0.)
Out[3]:
In [4]:
# Import TEAK NDVI Geotif
ndvi_filename = 'TEAK_NDVI.tif'
ndvi_dataset = gdal.Open(ndvi_filename)
#Display the dataset dimensions:
rows_ndvi = ndvi_dataset.RasterYSize; print('x:',rows_ndvi)
cols_ndvi = ndvi_dataset.RasterXSize; print('y:',cols_ndvi)
print('NDVI Raster Info:')
print('driver:',ndvi_dataset.GetDriver().LongName)
print('bands:',ndvi_dataset.RasterCount)
print('geotransform:',ndvi_dataset.GetGeoTransform())
print('projection:',ndvi_dataset.GetProjection())
ndvi_mapinfo = ndvi_dataset.GetGeoTransform()
xMin = ndvi_mapinfo[0]
yMax = ndvi_mapinfo[3]
xMax = xMin + ndvi_dataset.RasterXSize/ndvi_mapinfo[1] #divide by pixel width
yMin = yMax - ndvi_dataset.RasterYSize/ndvi_mapinfo[5] #divide by pixel height (note sign +/-)
ndvi_ext = (xMin,xMax,yMin,yMax)
print('raster extent:',aspectNS_ext)
#Convert CHM raster to NumPy array
ndvi_raster = ndvi_dataset.GetRasterBand(1)
noDataVal = ndvi_raster.GetNoDataValue()
scaleFactor = ndvi_raster.GetScale()
print('no data value:',noDataVal)
print('scale factor:',scaleFactor)
ndvi_stats = ndvi_raster.GetStatistics(True,True)
print('CHM Statistics: Minimum=%.2f, Maximum=%.2f, Mean=%.3f, StDev=%.3f' %
(ndvi_stats[0], ndvi_stats[1], ndvi_stats[2], ndvi_stats[3]))
ndvi_array = ndvi_dataset.GetRasterBand(1).ReadAsArray(0,0,cols_ndvi,rows_ndvi).astype(np.float)
print(ndvi_array)
In [5]:
# Plot NDVI array
# from matplotlib import colors
import copy
# ndvi_array = copy.copy(ndvi_array)
# ndvi_array = ndvi_array[::-1] #?? do raster arrays need to be inversed to look right as geotif?
fig, ax = plt.subplots()
plt.imshow(ndvi_array,extent=ndvi_ext)
plt.colorbar(); plt.set_cmap('RdYlGn');
plt.title('TEAK NDVI')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
In [6]:
#Plot a histogram of NDVI values with 20 bins
plt.hist(ndvi_array.flatten(),20)
Out[6]:
In [7]:
#Mask out pixels that are north facing:
import copy
import numpy.ma as ma
#first copy the ndvi array so we can further select a subset
ndvi_gtpt6 = copy.copy(ndvi_array)
ndvi_gtpt6[ndvi_array<0.6]=np.nan
print(ndvi_gtpt6)
plt.imshow(ndvi_gtpt6,extent=ndvi_ext)
plt.colorbar(); plt.set_cmap('RdYlGn');
plt.title('TEAK NDVI > 0.6')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
In [8]:
#Now include additional requirement that slope is North-facing (i.e. aspectNS_array = 1)
ndvi_gtpt6_Nslope = copy.copy(ndvi_gtpt6)
ndvi_gtpt6_Nslope[aspectNS_array!=1]=np.nan
plt.imshow(ndvi_gtpt6_Nslope,extent=ndvi_ext)
plt.colorbar(); plt.set_cmap('RdYlGn');
plt.title('TEAK, North Facing & NDVI > 0.6')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
In [12]:
# %load ../hyperspectral_hdf5/array2raster.py
"""
Array to Raster Function from https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html)
"""
import gdal, osr #ogr, os, osr
import numpy as np
def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,epsg):
cols = array.shape[1]
rows = array.shape[0]
originX = rasterOrigin[0]
originY = rasterOrigin[1]
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(epsg)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
In [14]:
epsg = 32611 #WGS 84, UTM Zone 11N
array2raster('TEAK_Nslope_ndvi_gtpt6.tif',(xMin,yMax),1,-1,ndvi_gtpt6_Nslope,epsg)