I'm performing NDVI calculation on a Planet Scope 4 band image as per Planet's documentation
The following block of code is what I wrote:
Extract band data from original image in working directory
import rasterio import numpy
image_file = "20170430_194027_0c82_3B_AnalyticMS"
with rasterio.open(image_file) as src: band_red = src.read(3)
with rasterio.open(image_file) as src: band_nir = src.read(4)
from xml.dom import minidom
xmldoc = minidom.parse("20170430_194027_0c82_3B_AnalyticMS_metadata") nodes = xmldoc.getElementsByTagName("ps:bandSpecificMetadata")
Extract TOA correction coefficients from metadata file in directory
TOA_coeffs = {} for node in nodes: bn = node.getElementsByTagName("ps:bandNumber")[0].firstChild.data if bn in ['1', '2', '3', '4']:
i = int(bn)
value = node.getElementsByTagName("ps:ReflectanceCoefficient")[0].firstChild.data
TOA_coeffs[1] = float(value)
Calculate NDVI and save file
band_red = band_red * TOA_coeffs[3] band_nir = band_nir * TOA_coeffs[4]
numpy.seterr(divide = 'ignore', invalid = 'ignore')
NDVI = (band_nir.astype(float) - band_red.astype(float))/(band_nir + band_red) numpy.nanmin(NDVI), numpy.nanmax(NDVI)
kwargs = src.meta kwargs.update(dtype=rasterio.float32,
count = 1)
with rasterio.open('ndvi.tif', 'W', **kwargs) as dst: dst.write_band(1, NDVI.astype(rasterio.float32))
Add symbology and plot color bar
import matplotlib.pyplot as plt import matplotlib.colors as colors
class MidpointNormalize(colors.Normalize): def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
self.midpoint = midpoint
colors.Normalize.__init__(self, vmin, vmax, clip)
def __call__(self, value, clip=None):
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
return numpy.ma.masked_array(numpy.interp(value, x, y), >numpy.isnan(value))
min = numpy.nanmin(NDVI) min = numpy.nanmax(NDVI) mid = 0.1
fig = plt.figure(figsize= (20,10)) ax = fig.add_subplot(111)
cmap = plt.cm.RdYlGn
cax = ax.imshow(NDVI, cmap=cmap, clim=(min,max),
>norm=MidpointNormalize(midpoint=mid, vmin=min, vmax=max))
ax.axis('off') ax.set_title('NDVI_test', fontsize= 18, fontweight='bold')
cbar = fig.colorbar(cax, orientation= 'horizontal', shrink=0.65)
fig.savefig("output/NDVI_test.png", dpi=200, bbox_inches='tight',
>pad_inches=0.7)
plt.show()
Plot histogram for NDVI pixel value distribution
fig2 = plt.figure(figsize=(10,10)) ax = fig2.add_subplot(111)
plt.title("NDVI Histogram", fontsize=18, fontweight='bold') plt.xlabel("NDVI values", fontsize=14) plt.ylabel("# pixels", fontsize=14)
x = NDVI[~numpy.isnan(NDVI)] numBins = 20 ax.hist(x,numBins,color='green',alpha=0.8)
fig2.savefig("output/ndvi-histogram.png", dpi=200, bbox_inches='tight', >pad_inches=0.7)
plt.show()
Alas, the execution of the script is cut short at the beginning of the code:
File "C:/Users/David/Desktop/ArcGIS files/Planet Labs/2017.6_Luis_Bedin_Bolivia/planet_order_58311/20170430_194027_0c82/TOA_correction_NDVI.py", line 8, in <module>
import rasterio
ModuleNotFoundError: No module named 'rasterio'
So I decide to install rasterio, that should solve the problem:
C:\Users\David\Desktop\ArcGIS files\Planet Labs\2017.6_Luis_Bedin_Bolivia\planet_order_58311\20170430_194027_0c82>pip install rasterio
Collecting rasterio
Using cached rasterio-0.36.0.tar.gz
Requirement already satisfied: affine in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: cligj in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: numpy in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: snuggs in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: click-plugins in c:\users\david\anaconda3\lib\site-packages (from rasterio)
What I interpret from this is that rasterio is already installed. How can this be if the Python console tells me there's no module named rasterio. The output from the console also says Microsoft Visual C++ is required. Upon further research I find this user's solution. I tried it but the console also tells me that rasterio is already installed:
(envpythonfs) C:\Users\David\Desktop\ArcGIS files\Planet Labs\2017.6_Luis_Bedin_Bolivia\planet_order_58311\20170430_194027_0c82>conda install rasterio gdal
Fetching package metadata .............
Solving package specifications: .
# All requested packages already installed.
# packages in environment at C:\Users\David\Anaconda3\envs\envpythonfs:
#
I'm creating the script using Spyder 3.1.2 with Python 3.6 on a Windows 10 64-bit machine.