Note
Click here to download the full example code
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
[[160 140 115]
[207 187 162]
[193 173 148]
...
[136 109 100]
[108 75 68]
[130 135 112]]
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>.
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.
Run the script
rM.run()
print(time.time()-t)
Out:
rasterMath... [........................................]0%
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%
Saved /tmp/sub.tif using function sub
0.3649485111236572
Plot result
#from osgeo import gdal
#from matplotlib import pyplot as plt
#src = gdal.Open('/tmp/sub.tif')
#plt.imshow(src.ReadAsArray())
Total running time of the script: ( 0 minutes 0.391 seconds)