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
[[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[-- -- --]
[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())
Out:
<matplotlib.image.AxesImage object at 0x7f7fdf6a4450>
Total running time of the script: ( 0 minutes 1.284 seconds)