Cómo reclasificar un ráster basado en cuantiles y descartando los no data values con Python

En el código siguiente los nodata values del array de un ráster se convierten en valores NaN (Not a Number) para aprovechar el método 'nanpercentile' de numpy que permite determinar estos estadísticos fácilmente excluyendo estos valores (algo que no se produce con 'percentile' y los nodata values).

Los percentiles determinados de esta manera se usan para reclasificar los valores del array según las condiciones establecidas en el código. Los nodata values del ráster original se recuperan nuevamente a partir de los NaN antes de grabar el ráster resultante con todos los parámetros del ráster original (geotransform, nodata values y proyección espacial).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64

import numpy as np from osgeo import gdal, gdal_array, osr dataset = gdal.Open("/home/zeito/pyqgis_data/test_raster_nodata.tif") band = dataset.GetRasterBand(1) nodata = band.GetNoDataValue() array = dataset.ReadAsArray() new_array = array nan_array = array nan_array[array == nodata] = np.nan percentile_80 = np.nanpercentile(nan_array, 80) percentile_60 = np.nanpercentile(nan_array, 60) percentile_40 = np.nanpercentile(nan_array, 40) percentile_20 = np.nanpercentile(nan_array, 20) percentile_0 = np.nanpercentile(nan_array, 0) for i, v in enumerate(new_array): for j, element in enumerate(v): if element <= percentile_0: new_array[i,j] = 1 if element > percentile_0 and element <= percentile_20: new_array[i,j] = 2 if element > percentile_20 and element <= percentile_40: new_array[i,j] = 3 if element > percentile_40 and element <= percentile_60: new_array[i,j] = 4 if element > percentile_60 and element <= percentile_80: new_array[i,j] = 5 if element > percentile_80: new_array[i,j] = 6 new_array[ new_array != new_array ] = nodata geotransform = dataset.GetGeoTransform() wkt = dataset.GetProjection() # Create gtif file driver = gdal.GetDriverByName("GTiff") output_file = "/home/zeito/pyqgis_data/test_raster_nodata_reclass.tif" dst_ds = driver.Create(output_file, band.XSize, band.YSize, 1, gdal.GDT_Int16) #writting output raster dst_ds.GetRasterBand(1).WriteArray( new_array ) #setting nodata value dst_ds.GetRasterBand(1).SetNoDataValue(nodata) #setting extension of output raster # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution dst_ds.SetGeoTransform(geotransform) # setting spatial reference of output raster srs = osr.SpatialReference() srs.ImportFromWkt(wkt) dst_ds.SetProjection( srs.ExportToWkt() ) #Close output raster dataset dataset = None dst_ds = None

Después de ejecutar el código anterior con el ráster siguiente; donde los valores más alla de la máscara del ráster (10x12) son nodata values establecidos como -999.

se obtiene el ráster de la imagen a continuación donde con el Value Tool plugin se corroboraron que sus values estaban acorde con la clasificación por cuantiles establecida en el código.