2024-05-17

PyQGISでNDVI

まだうまく動いてない
NIRとRのバンドが取り出せてない

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

### READ FROM RASTER

rasterPath = "/Users/shohei/Desktop/Drone/103FPLAN/New Folder With Items/_stack_raster.tif"
layer = iface.addRasterLayer(rasterPath,"my raster")
ext = layer.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
print(xmin,xmax,ymin,ymax)

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

#NDVI = (G-NIR)/(G+NIR)
#NDVI = (NIR-R)/(NIR+R)


### WRITE TO RASTER

src = rasterio.open(rasterPath)
fn = '/Users/shohei/Downloads/ndvi.tif'
driver = gdal.GetDriverByName('GTiff')
ds_write = driver.Create(fn, xsize=src.width, ysize=src.height, bands=1)

ds_write.GetRasterBand(1).WriteArray(NDVI)
#ds_write.GetRasterBand(2).WriteArray(band2)
#ds_write.GetRasterBand(3).WriteArray(band3)

geot = [xmin,(xmax-xmin)/float(src.width), 0, ymax, 0, (ymin-ymax)/float(src.height)]
ds_write.SetGeoTransform(geot)
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:32737")
ds_write.SetProjection(srs.ExportToWkt())
ds_write = None

rlayer = iface.addRasterLayer(fn, "NDVI")

RGB　この時点でおかしい

NDVI

2024-05-17

PyQGISでOpenCVを用いて画像処理

以下のgeotの計算方法
geot = [260270.3168900000164285,116.45/5000.0, 0, 9856611.2914200015366077, 0, -116.45/5000.0]

ラプラシアンフィルタの例

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

### READ FROM RASTER

rasterPath = "/Users/shohei/Downloads/Eunice_transparent_mosaic_group1_3_5.tif"
layer = iface.addRasterLayer(rasterPath,"my raster")

ds_read = gdal.Open(rasterPath)
band1 = np.array(ds_read.GetRasterBand(1).ReadAsArray())
band2 = np.array(ds_read.GetRasterBand(2).ReadAsArray())
band3 = np.array(ds_read.GetRasterBand(3).ReadAsArray())

### COMPUTER VISION MANIPULATION

band1 = cv2.Laplacian(band1,cv2.CV_64F)
band2 = cv2.Laplacian(band2,cv2.CV_64F)
band3 = cv2.Laplacian(band3,cv2.CV_64F)

### WRITE TO RASTER

src = rasterio.open(rasterPath)
fn = '/Users/shohei/Downloads/laplacian.tif'
driver = gdal.GetDriverByName('GTiff')
ds_write = driver.Create(fn, xsize=src.width, ysize=src.height, bands=3)

ds_write.GetRasterBand(1).WriteArray(band1)
ds_write.GetRasterBand(2).WriteArray(band2)
ds_write.GetRasterBand(3).WriteArray(band3)

geot = [260270.3168900000164285,116.45/5000.0, 0, 9856611.2914200015366077, 0, -116.45/5000.0]
ds_write.SetGeoTransform(geot)
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:32737")
ds_write.SetProjection(srs.ExportToWkt())
ds_write = None

rlayer = iface.addRasterLayer(fn, "laplacian")

ソーベルフィルタの例

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

### READ FROM RASTER

rasterPath = "/Users/shohei/Downloads/Eunice_transparent_mosaic_group1_3_5.tif"
layer = iface.addRasterLayer(rasterPath,"my raster")

ds_read = gdal.Open(rasterPath)
band1 = np.array(ds_read.GetRasterBand(1).ReadAsArray())
band2 = np.array(ds_read.GetRasterBand(2).ReadAsArray())
band3 = np.array(ds_read.GetRasterBand(3).ReadAsArray())

### COMPUTER VISION MANIPULATION

src_median = cv2.medianBlur(band1, 5)
sobel_x = cv2.Sobel(src_median, cv2.CV_32F, 1, 0) # X方向
sobel_y = cv2.Sobel(src_median, cv2.CV_32F, 0, 1) # Y方向
sobel_x = cv2.convertScaleAbs(sobel_x, alpha = 0.5)
sobel_y = cv2.convertScaleAbs(sobel_y, alpha = 0.5)
band1 = cv2.add(sobel_x, sobel_y)

src_median = cv2.medianBlur(band2, 5)
sobel_x = cv2.Sobel(src_median, cv2.CV_32F, 1, 0) # X方向
sobel_y = cv2.Sobel(src_median, cv2.CV_32F, 0, 1) # Y方向
sobel_x = cv2.convertScaleAbs(sobel_x, alpha = 0.5)
sobel_y = cv2.convertScaleAbs(sobel_y, alpha = 0.5)
band2 = cv2.add(sobel_x, sobel_y)

src_median = cv2.medianBlur(band3, 5)
sobel_x = cv2.Sobel(src_median, cv2.CV_32F, 1, 0) # X方向
sobel_y = cv2.Sobel(src_median, cv2.CV_32F, 0, 1) # Y方向
sobel_x = cv2.convertScaleAbs(sobel_x, alpha = 0.5)
sobel_y = cv2.convertScaleAbs(sobel_y, alpha = 0.5)
band3 = cv2.add(sobel_x, sobel_y)

### WRITE TO RASTER

src = rasterio.open(rasterPath)
fn = '/Users/shohei/Downloads/sobel.tif'
driver = gdal.GetDriverByName('GTiff')
ds_write = driver.Create(fn, xsize=src.width, ysize=src.height, bands=3)

ds_write.GetRasterBand(1).WriteArray(band1)
ds_write.GetRasterBand(2).WriteArray(band2)
ds_write.GetRasterBand(3).WriteArray(band3)

geot = [260270.3168900000164285,116.45/5000.0, 0, 9856611.2914200015366077, 0, -116.45/5000.0]
ds_write.SetGeoTransform(geot)
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:32737")
ds_write.SetProjection(srs.ExportToWkt())
ds_write = None

rlayer = iface.addRasterLayer(fn, "sobel")

2024-05-13

macOSのQGISでSemi Automatic Classification (SAC)プラグインのエラー

rs_calcが無いといわれるので、強制的にインポートした
/Users/shohei/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/SemiAutomaticClassificationPlugin/interface/band_calc_tab.py

 #try:
 #    from remotior_sensus.tools import band_calc as rs_calc
 #except Exception as error:
 #    str(error)
 #明示的にインポート
 from remotior_sensus.tools import band_calc as rs_calc

2024-05-12

macOSでGRASSが起動しない

ロケールの設定をGrass.shに追加したら起動した

/Applications/GRASS-8.3.app/Contents/MacOS/Grass.sh

  1 #! /bin/bash
  2 #############################################################################
  3 #
  4 # MODULE:     GRASS Initialization
  5 # AUTHOR(S):  Justin Hickey - Thailand - jhickey@hpcc.nectec.or.th
  6 #             William Kyngesburye - kyngchaos@kyngchaos.com
  7 #             Eric Hutton
  8 #             Michael Barton - michael.barton@asu.edu
  9 # PURPOSE:    The source file for this shell script is in
 10 #           macosx/app/grass.sh.in and is the grass startup script for
 11 #               the Mac OS X application build.
 12 # COPYRIGHT:    (C) 2000-2018 by the GRASS Development Team
 13 #
 14 #               This program is free software under the GNU General Public
 15 #           License (>=v2). Read the file COPYING that comes with GRASS
 16 #           for details.
 17 #
 18 #############################################################################
 19
## 以下を追加　
 20 export LC_ALL=en_US.UTF-8
 21 export LANG=en_US.UTF-8