Raszter réteg vizsgálata

Töltsük be a következő rétegeket:

  • Sentinel-2 műholdfelvételek csempe (RGB, MS), 2019

  • Települések

Töltsük be a Sentinel-2 csempe egy sávját, és vizsgáljuk meg a térbeli kiterjedését, jellemzőit:

import os
import processing

path = "d:/indikatrix-org/raszter-reteg-vizsgalata/" 
imgfolder = path + "sentinel2a/"

files = [imgfolder + name for name in os.listdir(imgfolder) if 
         os.path.isfile(os.path.join(imgfolder, name)) and os.path.join(
             imgfolder, name).endswith(".jp2") ]

band1 = files[1]
rlayer = QgsRasterLayer(files[1], 'sentinel2a-band1')
rlayer.extent().toString()
# get the raster type: 0 = GrayOrUndefined (single band), 1 = Palette (single band),
# 2 = Multiband
rlayer.rasterType()

Vágjuk ki a rasztert a Balatonalmádi járás területére:

import os
import processing

path = "d:/indikatrix-org/raszter-reteg-vizsgalata/" 
imgfolder = path + "sentinel2a/"

files = [imgfolder + name for name in os.listdir(imgfolder) if
         os.path.isfile(os.path.join(imgfolder, name)) and os.path.join(
             imgfolder, name).endswith(".jp2") ]

band1 = files[1]

filename = "balatonalmadi-jaras_utm"
masklayer_utm = QgsVectorLayer(path + filename + ".shp", filename, "ogr")
outputpath =  path + filename + band1[-8:-4] + ".tif"

parameters = {'INPUT': band1,
            'MASK': masklayer_utm,
            'NODATA': -9999,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': outputpath }
processing.run('gdal:cliprasterbymasklayer', parameters)

Alakítsuk át a kódot úgy, hogy az több sávon is végrehajtsa a kivágást:

import os
import processing

path = "d:/indikatrix-org/raszter-reteg-vizsgalata/" 
imgfolder = path + "sentinel2a/"

files = [imgfolder + name for name in os.listdir(imgfolder) if
         os.path.isfile(os.path.join(imgfolder, name)) and os.path.join(
             imgfolder, name).endswith(".jp2") ]

filename = "balatonalmadi-jaras_utm"
masklayer_utm = QgsVectorLayer(path + filename + ".shp", filename, "ogr")

bandlist = [ files[3], files[2], files[1]]
for band in bandlist:
    outputpath =  path + filename + band[-8:-4] + ".tif"
    parameters = {'INPUT': band,
            'MASK': masklayer_utm,
            'NODATA': -9999,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': outputpath }
    processing.run('gdal:cliprasterbymasklayer', parameters)

Füzzük össze a sávokat, készítsünk egy többsávos RGB réteget:

imgfolder = path
files = [imgfolder + name for name in os.listdir(imgfolder) if
         os.path.isfile(os.path.join(imgfolder, name)) and os.path.join(
             imgfolder, name).endswith(".tif") ]
bandlist = [ files[2], files[1], files[0]]

#processing.algorithmHelp("gdal:merge")
outputpath =  files[0][:-8] + "_rgb" + ".tif"
parameters = {'INPUT': bandlist,
            'PCT': False,
            'SEPARATE': True,
            'OUTPUT': outputpath }
processing.run('gdal:merge', parameters)

Készítsük el a járás településeihez tartozó színes műholdfelvétel kivágatokat:

filename = "balatonalmadi-telepulesek_utm"
masklayer_utm = QgsVectorLayer(path + filename + ".shp", filename, "ogr")
QgsProject.instance().addMapLayer(masklayer_utm, False)
rgb_layer = QgsRasterLayer(path + "balatonalmadi-jaras_utm__rgb.tif", 'jaras-rgb')

crs = masklayer_utm.crs()
idx = masklayer_utm.fields().lookupField('NAME') 
settlements = masklayer_utm.getFeatures()


for settlement in settlements:
    id = [settlement.id()]
    masklayer_utm.select(id)
    
    settlement_name = settlement.attributes()[idx]
    outputpath =  path + "clips/" + settlement_name + "_rgb" + ".tif"
    parameters = {
            'INPUT': rgb_layer,
            'MASK': QgsProcessingFeatureSourceDefinition(masklayer_utm.id(), True),
            'SOURCE_CRS': masklayer_utm.crs(),
            'TARGET_CRS': masklayer_utm.crs(),
            'NODATA': -9999,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': outputpath 
            }
    processing.run('gdal:cliprasterbymasklayer', parameters)    
    masklayer_utm.removeSelection()
    templayer = QgsRasterLayer(outputpath, settlement_name)
    QgsProject.instance().addMapLayer(templayer, True)

1. ábra: Balatonalmádi járás településeinek műholdfelvétel kivágatai.
1. ábra: Balatonalmádi járás településeinek műholdfelvétel kivágatai.

Megjegyzés: Feltűnik-e a hiba?

Bővítsük ki a kódot úgy, hogy mentsük el az egyes sávok hisztogramjaihoz szükséges adatokat:

    for i in range(1,4):
        geojson = path + "clips/" + settlement_name + '_' + str(i) + ".geojson"
        parameters = {'INPUT': templayer,
            'BAND': i, 
            'OUTPUT_TABLE': geojson }
        processing.run('native:rasterlayeruniquevaluesreport', parameters)

Illetve töltsük be az adatokat (JSON), és készítsünk hisztogramot, mentsük fájlba:

import matplotlib.pyplot as plt
import json
path = "d:/indikatrix-org/raszter-reteg-vizsgalata/" 
settlement_name = "Balatonkenese"

for i in range(1,4):
    settlement_json = path + "clips/" + settlement_name + '_' + str(i) + ".geojson"
    with open(settlement_json) as f:
        data = json.load(f)
    values = []
    counts = []
    for json_feature in data['features']:
        values.append(int(json_feature['properties']['value']))
        counts.append(int(json_feature['properties']['count']))
    #https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.hist.html
    plt.hist(values, bins = 100, weights=counts, histtype = 'step')
    #plt.hist(values, bins=np.arange(0,5050,50), weights=counts, histtype = 'step')
plt.show()
plt.savefig(path + "clips/" + "Balatonkenese-rgb-hist.png")

A kódrészlet beágyazását, továbbfejlesztését a kedves olvasóra bízom.

2. ábra: Balatonkenese RGB hisztogramja.
2. ábra: Balatonkenese RGB hisztogramja.

© Dr. Wirth Ervin