numpyRasterCalc.py

#########################
# NumPy Raster Calculator
# -----------------------
# Purpose: The purpose of this Python script is to use ArcGIS's arcpy package along with the Python library numpy to convert Landsat thermal bands from digital numbers to Top of Atmosphere (ToA) brightness-temperature values that provide pixel temperature in degrees Celsius.
# Purpose #2: The second purpose is that it allows some work to be done on rasters without the need for the Spatial Analyst extension
#
# Requirements: ESRI ArcGIS arcpy; Python 2.7.x; numpy; Landsat data
# ArcPy takes care of the dirty work of conversions and geospatial stuff and I have ArcGIS installed, hence it's use.
#########################

# Where code for the initial structure of this script came from: http://blog.remotesensing.io/2013/05/using-arcpy-for-raster-analysis/
#   This script uses the calculations found here (same as the USGS instructions, but in simpler form): http://geohackers.in/2013/08/using-data-from-the-landsat-8-tirs-instrument-to-estimate-surface-temperature/

# Maybe I'll also make another version of this script using Rasterio? I don't know... https://github.com/mapbox/rasterio

import arcpy
import numpy as np


infile = "D:\\GIS\\Landsat\\L8_S_20130715\\band10.tif"
in_arr = arcpy.RasterToNumPyArray(infile)#, nodata_to_value=9999)

##
#   Start Calculations
##

#Constants for the band and scene
rescaleFactorMult = 0.0003342   #band-specific multiplicative rescaling factor (Ml)
rescaleFactorAdd = 0.1          #band-specific additive rescaling factor (Al)
k1 = 774.89                     #band specific thermal conversion constant
k2 = 1321.08                    #band specific thermal conversion constant

print('Starting raster calculations')

#Digital Number to ToA Radiance Values
print('Calculating ToA spectral radiance values (converting digital number to radiance)')
toaSpectRad = (rescaleFactorMult * in_arr) + rescaleFactorAdd

#ToA Radiance to ToA Brightness-Tempearture in Kelvin
print('ToA Radiance to ToA Brightness-Tempearture in Kelvin')
toaBrightTemp = k2 / (np.log((k1/toaSpectRad + 1)))

#Kelvin to degrees Celsius
print('Converting Kelvin values to degrees Celsius')
out_arr = toaBrightTemp - 273.16

##
#   End Calculations
##

spatialref = arcpy.Describe(infile).spatialReference
cellsize1 = arcpy.Describe(infile).meanCellHeight
cellsize2 = arcpy.Describe(infile).meanCellWidth
extent = arcpy.Describe(infile).Extent
pnt = arcpy.Point(extent.XMin, extent.YMin)

#Save the raster
print('Converting the number array back into a raster')
out_ras = arcpy.NumPyArrayToRaster(out_arr,pnt,cellsize1,cellsize2)
print('Copying the raster to the directory')
arcpy.CopyRaster_management(out_ras,"D:\\GIS\\Landsat\\L8_S_20130715\\outputThermal.tif")
print("Defining the raster's projection (the same as the input)")
arcpy.DefineProjection_management("D:\\GIS\\Landsat\\L8_S_20130715\\outputThermal.tif", spatialref)