geologyidea.com
Viac

Extrakcia pixelového podielu tried jemnejšieho rastra v hrubšom

Extrakcia pixelového podielu tried jemnejšieho rastra v hrubšom


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.


Mám dva rastre, jeden má rozlíšenie 1 km. Druhý má priestorové rozlíšenie 30 m. Táto druhá v pixeloch predstavuje klasifikáciu. Chcem prekryť 1 km raster na klasifikačnom rastri a extrahovať pre každý 1 km pixel rozloženie tried zodpovedajúce 30 m klasifikácii.

Klasifikácia 30 m zahŕňa 7 tried, takže by som chcel, aby tento proces vytvoril 7 nových vrstiev na 1 km. Pre každú triedu jeden. Každý pixel 1 km bude mať teraz priradené zodpovedajúce% prítomnosti každej triedy prítomnej v klasifikácii 30 m. Ak by niekto vytvoril tabuľku, v ktorej je každý riadok pixelom, mal by tabuľku, ako je táto:

Ďalšou ťažkosťou je, že rastre sú dosť veľké. Naozaj by som to chcel dosiahnuť pomocou pythonu.

ÚPRAVA:

Pridal som podmnožiny svojich dvoch súborov údajov na tento odkaz (v rozlíšení 30 ma 1000 m), pre prípad, že by niekto chcel navrhnúť dôkaz koncepcie; bolo by úžasné.


Krok 1

Vytvorte bitové rastre pre každú z jedinečných tried. Môže to byť 1-pásmový raster pre každú triedu alebo jeden raster s pásmom pre každú triedu (napr. GeoTIFF). Ak používate GTiff, môžete použiť možnosť vytvoreniaNBITS = 1šetriť priestor. Môžete tiež zvážiť použitie twobitových rastrov na ukladanie logiky s tromi hodnotami, kde tretí (napr. 2) je NODATA, čo by vyžadovaloNBITS = 2.

Pre každé pásmo / triedu budú pixely 0, ak táto trieda neexistuje, alebo 1, ak táto trieda existuje.

import numpy ako np z osgeo import gdal gdal.UseExceptions () # Získanie údajov z rastra s klasifikáciami ds = gdal.Open ('cropped_30m.tif') band = ds.GetRasterBand (1) class_ar = band.ReadAsArray () gt = ds .GetGeoTransform () pj = ds.GetProjection () ds = pásmo = žiadne # zavrieť # Definujte rastrové hodnoty pre každú triedu, aby sa vzťahovali ku každému pásmu class_ids = (np.arange (31) + 1) .tolist () # Make nový bitový raster drv = gdal.GetDriverByName ('GTiff') ds = drv.Create ('bit_raster.tif', class_ar.shape [1], class_ar.shape [0], len (class_ids), gdal.GDT_Byte, [ 'NBITS = 1', 'COMPRESS = LZW', 'INTERLEAVE = BAND']) ds.SetGeoTransform (gt) ds.SetProjection (pj) pre bidx v rozsahu (ds.RasterCount): band = ds.GetRasterBand (bidx + 1 ) # create boolean selection = (class_ar == class_ids [bidx]) band.WriteArray (selection.astype ('B')) ds = band = None # uložiť, zavrieť

Napr. pre červené = pásmo 10, zelené = pásmo 13 a modré = pásmo 27 je možné ho vizualizovať:

Čierna označuje, že ide o inú triedu / pásmo.

Krok 2

Použite gdalwarp na získanie priemer 0 a 1 hodnôt z hrubšej mriežky vzorky. Táto metóda vyžaduje GDAL> = 1.10.0. Skutočne jednoduchý spôsob, ako to dosiahnuť, je zadať cieľové rozlíšenie na príkazovom riadku, napr. 1 km:

$ gdalwarp -ot Float32 -tr 1000 1000 -r priemerný bit_raster.tif priemerný_1km.tif

Výsledok bude mať toľko pásiem ako tried a každé pásmo bude popisovať zlomok tejto triedy pre každý pixel v rozmedzí od 0 do 1. Ak dávate prednosť percentám, vynásobte ich číslom 100 (napr. Pozri gdal_calc.py).

Existuje tiež rozhranie Python pre gdalwarp (napr. Táto odpoveď), ktoré môže používať šablónový raster pre tvar a kde by sa mala použiť NODATA. Napríklad:

# Otvorte raster z kroku 1 src_ds = gdal.Open ('bit_raster.tif') # Otvorte šablónu alebo skopírujte pole pre rozmery a masku NODATA cpy_ds = gdal.Open ('cropped_1000m.tif') band = cpy_ds.GetRasterBand (1 ) cpy_mask = (band.ReadAsArray () == band.GetNoDataValue ()) # Výsledný raster, s rovnakým rozlíšením a pozíciou ako rastrový kopírovací dst_ds = drv.Create ('average2_1km.tif', cpy_ds.RasterXSize, cpy_ds.RasterYSize, len (class_ids), gdal.GDT_Float32, ['INTERLEAVE = BAND']) dst_ds.SetGeoTransform (cpy_ds.GetGeoTransform ()) dst_ds.SetProjection (cpy_ds.GetProjection ()) # Urobte to isté ako priemerná hodnota gdalwarp -r; dokončenie gdal.ReprojectImage (src_ds, dst_ds, None, None, gdal.GRA_Average) môže chvíľu trvať # Preveďte všetky zlomky na percentá a použite rovnakú masku NODATA z rastra kopírovania NODATA = -9999 pre bidx v rozsahu (dst_ds .RasterCount): band = dst_ds.GetRasterBand (bidx + 1) ar = band.ReadAsArray () * 100,0 ar [cpy_mask] = NODATA band.WriteArray (ar) band.SetNoDataValue (NODATA) # Uložiť a zavrieť všetky rastre src_ds = cpy_ds = dst_ds = band = žiadny

Napr. rovnaké tri pásma (triedy) je možné vizualizovať ako RGB:

Alebo ho stiahnite do QGIS a pomocou nástroja Identifikovať funkcie skontrolujte zlomok pre každé pásmo / klasifikáciu na jednom mieste:

Tieto čísla dosahujú až 100%, čo je dobrá kontrola.


Pozri si video: Výpočet čisté mzdy u DPP a DPČ