El blog de José Guerrero: Extraer el punto máximo de un ráster para una rejilla de 4x4 con los métodos de GDAL Python en QGIS 3

Esto puede realizarse con el método 'ReadAsArray' del módulo Python de GDAL empleando una ventana de 4x4 y almacenando cada punto donde ocurren los máximos en una memory layer. Para probar esta aproximación se creó un ráster de resolución de 0.5x0.5, con una amplitud (filas, columnas) de 20x20 y valores entre 1 y 1000 para evitar valores repetidos. El código desarrollado se expone a continuación y contiene varios print con propósito de corroboración.

from osgeo import gdal, osr
import numpy as np

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/test_0_5_0_5_clip.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

transform = dataset.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

cols = dataset.RasterXSize
rows = dataset.RasterYSize

nb = 4

idxs = []

for i in range(0, rows, nb):
    for j in range(0, cols, nb):
        data = band.ReadAsArray(j, i, nb, nb)
        print (data)
        max = np.max(data)
        print ("max: ", max)
        print ("indices maximum original matrix")

        try:
            for m, item in enumerate(data):
                for n, element in enumerate(item):
                    if element == max:
                        idxs.append([m + i, n + j])
                        print(m + i, n + j)

        except TypeError:
            print ("there is no matrix")

x0 = xOrigin + pixelWidth/2
y0 = yOrigin - pixelHeight/2

points = []

for item in idxs:
    x = x0 + item[1]*pixelWidth
    y = y0 - item[0]* pixelHeight
    points.append(QgsPointXY(x,y))

wkt = dataset.GetProjection()

crs = QgsCoordinateReferenceSystem()
crs.createFromWkt(wkt)

epsg = crs.postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPointXY(points[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

Al ejecutar el código anterior, se imprimen en la Python Console de QGIS 3 las sub matrices para las diferentes ventanas de 4x4 que permiten verificar, con el plugin Value Tool de QGIS, que los puntos generados en la memory layer corresponden a los máximos en cada ventana.

[[557 231 752 950]
 [937   3 282 817]
 [171 548 846 977]
 [929 788 445  64]]
max:  977
indices maximum original matrix
2 3
[[240 503   7 417]
 [ 84 691 877   7]
 [762 163 243 556]
 [644 822 484 860]]
max:  877
indices maximum original matrix
1 6
[[ 38 124 342 813]
 [435 548 524  62]
 [286 101 439 690]
 [598  64   6 589]]
max:  813
indices maximum original matrix
0 11
[[773 594 866 548]
 [427 921 883 621]
 [730 262 270 757]
 [464 158 197 969]]
max:  969
indices maximum original matrix
3 15
[[786 997 865 821]
 [166 118 365 254]
 [483 398 212 419]
 [116 349 555 835]]
max:  997
indices maximum original matrix
0 17
[[640 608 836 124]
 [583 300 605 725]
 [888 927 848 377]
 [866 858 101 396]]
max:  927
indices maximum original matrix
6 1
[[161 301 214 915]
 [760 707 844 739]
 [325 645 883 899]
 [120 307 909 711]]
max:  915
indices maximum original matrix
4 7
[[376 232 369 914]
 [454 300 472 453]
 [970 350 651 959]
 [920 102 706 305]]
max:  970
indices maximum original matrix
6 8
[[ 83  45 620 414]
 [553 493 856 315]
 [572 532 872 384]
 [162 176  91 839]]
max:  872
indices maximum original matrix
6 14
[[222 100   6 767]
 [626 286 227 360]
 [ 58 138 392 251]
 [295 870 719 885]]
max:  885
indices maximum original matrix
7 19
[[168 751 925 455]
 [950  27 809 331]
 [290 800 581 509]
 [ 92 400 144 934]]
max:  950
indices maximum original matrix
9 0
[[891 305 267 809]
 [ 19 457 455 554]
 [303 727 537 865]
 [ 66 253  40 740]]
max:  891
indices maximum original matrix
8 4
[[ 67 831 329 703]
 [804 990 563 263]
 [384 761 801 629]
 [849 441 620 479]]
max:  990
indices maximum original matrix
9 9
[[980 912 581 658]
 [375 879 202 317]
 [ 11 691 138 782]
 [548 538 145 217]]
max:  980
indices maximum original matrix
8 12
[[512 486 145 897]
 [ 79 549   7 618]
 [401 378 826 577]
 [944  35 716 875]]
max:  944
indices maximum original matrix
11 16
[[334 815 132  13]
 [103 767 803 610]
 [532 493 457 809]
 [420  11 311 198]]
max:  815
indices maximum original matrix
12 1
[[201 471 685 237]
 [670   8 707 373]
 [  8 580 690 447]
 [630 869   1   7]]
max:  869
indices maximum original matrix
15 5
[[108 490 358 911]
 [669 285 108 640]
 [148 843 546 293]
 [701 563  79 870]]
max:  911
indices maximum original matrix
12 11
[[339 768 761 774]
 [423  30 326 145]
 [ 22  32 270 764]
 [832 512 672 950]]
max:  950
indices maximum original matrix
15 15
[[416 194 733 435]
 [937 868 344 236]
 [343 852 713 597]
 [163 536 378 610]]
max:  937
indices maximum original matrix
13 16
[[ 64 780 655 329]
 [710 583  48 390]
 [256  73 133  96]
 [446 139 587 703]]
max:  780
indices maximum original matrix
16 1
[[522 609 373 742]
 [117 780 443   2]
 [151 755 598  50]
 [891 365 398 589]]
max:  891
indices maximum original matrix
19 4
[[986 237 501 787]
 [237 586 376 773]
 [187 597 884 789]
 [944 770 441 730]]
max:  986
indices maximum original matrix
16 8
[[915 868  32 780]
 [156 528 133  86]
 [707 945 425 766]
 [190   5 142 960]]
max:  960
indices maximum original matrix
19 15
[[942 404 576 803]
 [282 122 496 181]
 [892 781 924  54]
 [ 30 149 168 742]]
max:  942
indices maximum original matrix
16 16

En la imagen siguiente se pueden observar los puntos en el centro de cada píxel donde ocurren los máximos (también corroborado con el plugin Value Tool de QGIS). Se añadió una rejilla de 4x4 (2m x 2m) para facilitar la visualización.

El blog de José Guerrero

viernes, 22 de noviembre de 2019

Extraer el punto máximo de un ráster para una rejilla de 4x4 con los métodos de GDAL Python en QGIS 3

No hay comentarios: