Basics to use rasterMath

Compute substract and addition between two raster bands.

Import librairies

from museotoolbox.processing import RasterMath
from museotoolbox import datasets
import numpy as np

Load HistoricalMap dataset

raster,vector = datasets.load_historical_data()

Initialize rasterMath with raster

rM = RasterMath(raster)

print(rM.get_random_block())

Out:

Total number of blocks : 15
[[-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [198 178 145]
 [208 188 155]
 [208 188 155]
 [214 194 161]
 [211 191 158]
 [209 189 156]
 [215 195 162]
 [216 196 163]
 [205 185 152]
 [220 200 167]
 [196 176 143]
 [196 176 143]
 [207 187 154]
 [148 128 95]
 [99 79 46]
 [95 75 42]
 [195 175 140]
 [214 192 155]
 [217 194 160]
 [218 195 161]
 [218 195 161]
 [215 192 160]
 [210 187 156]
 [207 183 155]
 [207 183 155]
 [176 152 124]
 [195 172 141]
 [204 181 149]
 [206 183 149]
 [206 183 149]
 [214 192 155]
 [213 191 154]
 [204 182 145]
 [204 184 151]
 [200 177 145]
 [211 187 151]
 [210 184 147]
 [213 185 146]
 [213 185 146]
 [195 171 137]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]
 [-- -- --]]

Let’s suppose you want compute the difference between blue and green band. I suggest you to define type in numpy array to save space while creating the raster!

X = rM.get_random_block()

sub = lambda X : np.array((X[:,0]-X[:,1])).astype(np.int16)


rM.add_function(sub,out_image='/tmp/sub_lambda.tif')

Out:

Using datatype from numpy table : int16.
Detected 1 band for function <lambda>.
No data is set to : -32768.

Use a python function to use arguments

def sub(X,band1=0,band2=1):
    outX = np.array((X[:,band1]-X[:,band2])).astype(np.int16)
    return outX

We can add keyword argument in the addFunction. This function is going to substract band2 from band 1

import time
t=time.time()
rM = RasterMath(raster)
rM.add_function(sub,out_image='/tmp/sub.tif',band1=1,band2=0,compress='high')

Out:

Total number of blocks : 15
Using datatype from numpy table : int16.
Detected 1 band for function sub.
No data is set to : -32768.

Run the script

rM.run()
print(time.time()-t)

Out:

Batch processing (15 blocks using 5Mo of ram)

rasterMath... [##......................................]6%

rasterMath... [#####...................................]13%

rasterMath... [########................................]20%

rasterMath... [##########..............................]26%

rasterMath... [#############...........................]33%

rasterMath... [################........................]40%

rasterMath... [##################......................]46%

rasterMath... [#####################...................]53%

rasterMath... [########################................]60%

rasterMath... [##########################..............]66%

rasterMath... [#############################...........]73%

rasterMath... [################################........]80%

rasterMath... [##################################......]86%

rasterMath... [#####################################...]93%

rasterMath... [########################################]100%
1.1145141124725342

Plot result

from osgeo import gdal
from matplotlib import pyplot as plt
src = gdal.Open('/tmp/sub.tif')
plt.imshow(src.ReadAsArray())
rasterMath

Out:

<matplotlib.image.AxesImage object at 0x7f7fdf6a4450>

Total running time of the script: ( 0 minutes 1.284 seconds)

Gallery generated by Sphinx-Gallery