10. Querying Rasters¶

Below are sample scripts that can be used to query rasters

Querying Rasters Using an ArcPy Map Algebra Expression

# Find all values in an elevation raster that greater than 800 (meters).

import arcpy
from arcpy.sa import *
arcpy.env.overwriteOutput = True

# Specify the input raster
inRaster = "C:/workspace/dem"

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

# Make a map algebra expression and save the resulting raster
outRaster = Raster(inRaster) > 800
outRaster.save("C:/workspace/query800c")

# Check in the Spatial Analyst extension now that you're done
arcpy.CheckInExtension("Spatial")
print "Query successfully executed.

Querying a Raster using Arcpy with the “Con” Keyword

‘Con’ is abbreviation for condition and allows for the execution of conditional statements. The basic syntax when using “Con” is:

Con(in_conditonal_raster, in_true_raster_or_constant, {in_false_raster_or_constant}, {where_clause}).

Example (a)

import arcpy
from arcpy import env
from arcpy.sa import *

env.overwriteOutput = True

# Check out ArcGIS Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

env.workspace = "C:/workspace"

outElev1 = Con(dem > 700, 1, 0)
outElev1.save("C:/workspace/outElev1")

#or

# Execute Con using a map algebra expression instead of a where clause
outCon2 = Con(Raster("dem") > 700, "dem")

outCon2.save("C:/workspace/outCon2")

Example (b)

import arcpy
from arcpy import env
from arcpy.sa import *

env.overwriteOutput = True

# Check out ArcGIS Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

env.workspace = "C:/workspace"

outCon = Con("dem", "dem", "", "VALUE > 700")
outCon.save("C:/workspace/outcon.img")

Using the Con tool with 2 different Rasters

import arcpy
from arcpy import env
from arcpy.sa import *

env.overwriteOutput = True
env.workspace = "c:/workspace"

# Check out ArcGIS Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

inSlope = Slope("dem", "DEGREE", 0.3043)
inSlope.save("C:/workspace/inslope01")

# Execute Aspect
inAspect = Aspect("dem")

# Save the output
inAspect.save("C:/workspace/inaspect02")


outCon = Con(((inSlope >= 0) & (inSlope <= 10)) & ((inAspect >= 110) & (inAspect <= 160)), 1, "")

outCon.save("C:/workspace/con56")

print ("Query Executed Successfully")

Querying a Raster Using ExtractByAttributes

import arcpy
from arcpy import env
from arcpy.sa import *

env.overwriteOutput = True
env.workspace = "c:/workspace"

# Check out ArcGIS Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")


attExtract = ExtractByAttributes("inslope01", "VALUE > 2 & VALUE < 10")
attExtract.save("C:/workspace/extract01")

print "Query Executed Successfully"

Querying Rasters with Non- ESRI Libraries

#a. Find all places with elevation greater than 2500 ft (Method 1)

import gdal
import matplotlib.pyplot as plt

#Open raster and read number of rows, columns, bands
dataset = gdal.Open("C:/Users/Hugh/Desktop/N47E010.hgt")
band = dataset.GetRasterBand(1)

elevation = band.ReadAsArray(0,0,cols,rows)

a2500 = elevation > 2500

plt.imshow(a2500)
plt.show()

#b. Find all places with elevation greater than 2500 ft (Method 2)

import gdal
import matplotlib.pyplot as plt

#Open raster and read number of rows, columns, bands
dataset = gdal.Open("C:/Users/Hugh/Desktop/N47E010.hgt")
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

elevation = band.ReadAsArray(0,0,cols,rows)

outData = np.zeros((rows, cols))

for y in range(rows):
     for x in range(cols):
          if elevation[y, x] > 2500:
              outData[y, x] = 1
          else:
              outData[y, x] = 0

plt.imshow(outData)
plt.show()

10.1. Map Algebra¶

Map algebra involves performing mathematical operations on rasters, as shown in the example below. In ArcGIS, such operations can be performed manually using the Raster Calculator. Map algebra allows us to plug rasters into equations to perform modeling in a GIS environment.

1. Simple Map Algebra Operations

import arcpy
from arcpy import env
from arcpy.sa import *
try:
# Set environment settings
env.workspace = "C:/workspace"
env.overwriteOutput = True

# Check out ArcGIS Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

outRas2 = Raster("dem") * 8
outRas3 = outRas2 + Raster("dem") * 10
print "Calculation successfully completed"


except Exception as e:
     print e.message

2. Implementing Algebraic Equations

The equation below calculates topographic wetness index for an area. You can read about it here - You can read about the index here - http://soilandwater.bee.cornell.edu/tools/TI.pdf.

Links to an external site. The script below shows how to deal with equations. In ArcPy, slope maps in degrees must be converted to radians.

Topographic Wetness Index

import arcpy
from arcpy import env
from arcpy.sa import *
import math


# Set environment settings
env.workspace = "C:\workspace"
env.overwriteOutput = True

# Set local variables
inRaster = "slope1"
flowaccum1 = "flowaccum"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute Ln converting. All degrees are converted to degrees to radians
outLn = Ln(((Raster(flowaccum1) + 1)*10) / Tan((Raster(inRaster) * (math.pi / 180.0))))

# Save the output
outLn.save("C:\workspace\outln6")

print "Computation executed successfully"

3. Calculating NDVI with GDAL

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors


band4 = gdal.Open("/Users/semple/Desktop/LT05/LT05_L1TP_047027_20101006_20200823_02_T1_B4.TIF")
band3 = gdal.Open("/Users/semple/Desktop/LT05/LT05_L1TP_047027_20101006_20200823_02_T1_B3.TIF")


b3 = band3.GetRasterBand(1)
b4 = band4.GetRasterBand(1)


fig, ax = plt.subplots(figsize= (8,8))

b3_arr = b3.ReadAsArray().astype(np.float32)
b4_arr = b4.ReadAsArray().astype(np.float32)

#Allow division by zero
np.seterr(divide='ignore', invalid='ignore')

ndvi = (b4_arr - b3_arr) / (b4_arr + b3_arr)

plt.imshow(ndvi)
plt.colorbar(shrink=0.50)

plt.show()

ndvi.png
numply.seterror(). - https://numpy.org/doc/stable/reference/generated/numpy.seterr.html

4. Calculating NDVI with Rasterio

import numpy as np
import rasterio
from rasterio.plot import show

blueband  = rasterio.open("C:/Users/Hugh/Desktop/Landsat/Landsat/LT05_L1TP_019031_20111106_20160830_01_T1_B4.TIF")
greenband = rasterio.open("C:/Users/Hugh/Desktop/Landsat/Landsat/LT05_L1TP_019031_20111106_20160830_01_T1_B3.TIF")
mirband    = rasterio.open("C:/Users/Hugh/Desktop/Landsat/Landsat/LT05_L1TP_019031_20111106_20160830_01_T1_B2.TIF")

green = greenband.read(1).astype(float)
mir = mirband.read(1).astype(float)

np.seterr(divide='ignore', invalid='ignore')  # Allow division by zero

mndwi = np.empty(greenband.shape, dtype=rasterio.float32)  # Create empty matrix

check = np.logical_or(mir > 0.0, green > 0.0)  # Create check raster with True/False values

mndwi = np.where(check, (green - mir) / (green + mir), -999)  # Calculate MNDWI

water = np.where(mndwi > 0, 1, 0) # Set values above 0 as water and otherwise leave it at 0

show(water, cmap='plasma')

Readings

https://colekrehbiel.wordpress.com/2017/10/18/calculate-ndvi-from-sentinel-2a-data/

https://github.com/ckrehbiel/website/blob/master/Sentinel_NDVI.py

https://stackoverflow.com/questions/16411468/calculate-ndvi-using-python

Python GIS Tutorials

Navigation

Related Topics

10. Querying Rasters¶

10.1. Map Algebra¶