Saturday, 2 May 2015

numpy - How to filter no data value with GDAL?


I have a Tiff file that has no data value = -3.40282347e+38. I want to filter no data and get raw values with gdal. I tried with below code:


import numpy as np

from osgeo import gdal, gdal_array

dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()

lst = []
for v in array:
if v < 0:
lst.append(v)


print lst

but I received this error:


Traceback (most recent call last):
File "", line 1, in
File "/tmp/tmp3B0d5H.py", line 10, in
if v < 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

How can I do it?




Answer



If you want to filter no data and get raw values you need following code:


import numpy as np
from osgeo import gdal, gdal_array

dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()

lst = [ element for v in array for element in v if element > 0 ]


print lst

I tried it out with a raster with no data values equal to -999 and it works adequately.


Editing Note:


Based in your commentary, you don't need to filter these values. It's necessary to change no data values to NaN values for calculating directly np.nanpercentiles. Following code calculates this (percentile 50) for my raster with -999 no data values.


import numpy as np
from osgeo import gdal, gdal_array

dataset = gdal.Open("/home/zeito/pyqgis_data/test_raster_nodata.tif")
array = dataset.ReadAsArray()


nan_array = array

nan_array[array == -999] = np.nan

print np.nanpercentile(nan_array, 50)

Result is 6.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...