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
[[140 136 91]
 [132 132 82]
 [141 144 91]
 [140 145 89]
 [140 148 89]
 [140 148 89]
 [144 152 95]
 [149 157 100]
 [155 163 106]
 [158 167 110]
 [159 156 87]
 [161 158 89]
 [163 160 93]
 [166 163 96]
 [166 163 96]
 [168 163 97]
 [169 164 100]
 [169 164 100]
 [168 163 99]
 [164 159 95]
 [164 156 93]
 [165 157 94]
 [162 154 91]
 [162 154 91]
 [158 150 87]
 [160 152 89]
 [166 158 93]
 [161 156 88]
 [173 165 100]
 [169 159 97]
 [169 160 103]
 [176 168 119]
 [176 168 119]
 [149 145 107]
 [123 126 95]
 [80 88 64]
 [89 102 84]
 [128 145 135]
 [139 156 146]
 [135 151 140]
 [161 173 159]
 [161 173 159]
 [164 173 154]
 [156 160 137]
 [190 190 162]
 [171 170 142]
 [176 170 144]
 [174 163 141]
 [173 155 135]
 [177 154 136]
 [217 180 153]
 [143 102 74]
 [187 136 107]
 [182 128 100]
 [210 163 133]
 [205 171 143]
 [208 183 153]
 [190 163 134]
 [193 166 137]
 [193 166 137]
 [200 173 144]
 [193 166 137]
 [188 163 133]
 [191 168 137]
 [182 159 128]
 [217 194 163]
 [201 180 149]
 [200 179 148]
 [200 179 148]
 [216 195 164]
 [189 168 137]
 [206 186 153]
 [210 187 155]
 [174 150 116]
 [193 169 135]
 [200 176 142]
 [203 179 145]
 [203 179 145]
 [200 176 142]
 [186 162 128]
 [192 168 134]
 [188 164 130]
 [211 187 153]
 [198 172 139]
 [203 177 144]
 [178 152 119]
 [178 152 119]
 [183 152 123]
 [193 153 128]
 [203 161 137]
 [211 166 143]
 [216 171 148]
 [216 171 148]
 [194 152 128]
 [174 129 106]
 [185 143 119]
 [133 88 65]
 [131 89 65]
 [139 94 71]
 [145 103 79]]

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.0585341453552246

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 0x7f3dd0432450>