PyQGISでNDVI

from osgeo import _gdal as gdal
from osgeo import osr
import numpy as np
import rasterio
import cv2
import tempfile

### READ FROM RASTER
layer = iface.layerTreeView().selectedLayers()[0]
rasterPath = layer.source()
ext = layer.extent()
width = layer.width()
height= layer.height()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
crs_string = layer.crs().authid()
geot = [xmin,(xmax-xmin)/float(width), 0, ymax, 0, (ymin-ymax)/float(height)]

ds_read = gdal.Open(rasterPath)
B = np.array(ds_read.GetRasterBand(1).ReadAsArray())
G = np.array(ds_read.GetRasterBand(2).ReadAsArray())
R = np.array(ds_read.GetRasterBand(3).ReadAsArray())
NIR = np.array(ds_read.GetRasterBand(4).ReadAsArray())


### COMPUTER VISION MANIPULATION
NDVI = (NIR-R)/(NIR+R) * 255/2.0
#NDVI = (NIR-R)/(NIR+R)*255/1.0 #almost the same as above result

### WRITE TO RASTER

OUTPUT_LAYER_NAME = "NDVI"
filename = OUTPUT_LAYER_NAME+".tif"
filepath = os.path.join(tempfile.mkdtemp(), filename)
tmpfile_read = open(filepath, "w")
tmpfile_read.close()
fn = filepath

driver = gdal.GetDriverByName('GTiff')
ds_write = driver.Create(fn, xsize=width, ysize=height, bands=1)
ds_write.GetRasterBand(1).WriteArray(NDVI)
#ds_write.GetRasterBand(2).WriteArray(band2)
#ds_write.GetRasterBand(3).WriteArray(band3)

ds_write.SetGeoTransform(geot)
srs = osr.SpatialReference()
srs.SetFromUserInput(crs_string)
ds_write.SetProjection(srs.ExportToWkt())
ds_write = None

rlayer = iface.addRasterLayer(fn, OUTPUT_LAYER_NAME)

RGB

NDVI (Raster calculator)

NDVI (PyQGIS)

ほとんど結果は変わらないので動いている