7

I'm trying to identify null/0 values in each of my raster bands. How can I loop through all the values of each band? I've tried this code, straight from the gdal/python cookbook website but this just gives stats for each band as a whole and for some reason it actually doesn't print a minimum value of 0 even though my raster does have 0 values, it shows a minimum value of 1.

I want to print every single cell value. Is there a way to do this? 

from osgeo import gdal import sys src_ds = gdal.Open( "INPUT.tif" ) print "[ RASTER BAND COUNT ]: ", src_ds.RasterCount for band in range( src_ds.RasterCount ): band += 1 print "[ GETTING BAND ]: ", band srcband = src_ds.GetRasterBand(band) stats = srcband.GetStatistics( True, True ) print "[ STATS ] = Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( \ stats[0], stats[1], stats[2], stats[3] )
Improve this question

4 Answers 4

Reset to default
7

You may read it as array, using numpy:

from osgeo import gdal import sys import numpy as np src_ds = gdal.Open( "INPUT.tif" ) print "[ RASTER BAND COUNT ]: ", src_ds.RasterCount for band in range( src_ds.RasterCount ): band += 1 print "[ GETTING BAND ]: ", band srcband = src_ds.GetRasterBand(band) stats = srcband.GetStatistics( True, True ) print "[ STATS ] = Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( \ stats[0], stats[1], stats[2], stats[3] ) rast_array = np.array(src_ds.GetRasterBand(1).ReadAsArray()) print rast_array

Using a sample raster, the above code will return:

[ RASTER BAND COUNT ]: 1 [ GETTING BAND ]: 1 [ STATS ] = Minimum=1683.000, Maximum=1900.000, Mean=1820.854, StdDev=59.329 [[1900 1900 1898 1895 1892 1887 1879 1871 1863 1852 1845 1837 1824 1802 1743 1725 1713 1705 1699 1693 1687 1683] [1897 1896 1894 1892 1890 1884 1877 1869 1862 1854 1847 1838 1820 1800 1745 1729 1719 1712 1706 1701 1696 1695] [1892 1891 1890 1888 1885 1881 1875 1868 1861 1855 1849 1837 1817 1794 1747 1732 1725 1720 1714 1710 1707 1706] [1887 1885 1884 1882 1880 1878 1873 1867 1860 1855 1849 1833 1815 1789 1749 1738 1732 1728 1723 1720 1718 1715] [1882 1880 1878 1876 1875 1873 1871 1866 1861 1855 1849 1832 1817 1795 1756 1744 1740 1737 1733 1730 1728 1725] [1880 1877 1874 1873 1870 1868 1867 1865 1860 1855 1850 1841 1834 1817 1795 1769 1749 1746 1743 1740 1736 1731] [1880 1876 1873 1870 1869 1866 1863 1862 1859 1856 1852 1847 1843 1841 1824 1812 1802 1775 1758 1747 1740 1733] [1879 1876 1873 1870 1869 1866 1863 1860 1858 1855 1852 1850 1847 1843 1831 1819 1803 1782 1763 1747 1738 1730] [1879 1877 1874 1872 1869 1866 1864 1861 1858 1855 1852 1850 1850 1848 1836 1816 1794 1775 1754 1744 1736 1728] [1880 1877 1875 1872 1869 1867 1864 1862 1858 1854 1850 1848 1850 1850 1840 1806 1786 1767 1749 1742 1734 1726] [1881 1879 1876 1873 1870 1866 1864 1861 1857 1851 1843 1840 1841 1850 1827 1797 1782 1769 1752 1742 1733 1723] [1882 1879 1876 1873 1870 1867 1864 1861 1855 1848 1839 1835 1833 1836 1810 1794 1783 1771 1758 1747 1737 1729] [1882 1880 1876 1873 1869 1866 1862 1858 1854 1849 1838 1833 1826 1814 1800 1792 1782 1773 1762 1752 1742 1733] [1881 1878 1874 1870 1867 1863 1860 1856 1853 1849 1840 1835 1821 1813 1798 1790 1783 1774 1766 1757 1748 1738]]

(If you want to print each value separately, it could be easy to edit the code).

Improve this answer
2
  • Thanks, this is exactly what I wanted. I've never worked with numpy. So it's making a bunch of arrays made up of individual cell values? So the first number 1900, that is the value of a singe cell correct? How does it determine number of values in each array? Commented Dec 8, 2016 at 17:26
  • Yes. The first array store the cell values starting from the up-left of the raster, and so on. Commented Dec 8, 2016 at 17:39
4

The problem has been resolved in GDAL does not ignore NoData value

 f = gdal.Open("a.tif") bands = f.RasterCount print bands 3 for j in range(bands): band = f.GetRasterBand(j+1) stats = band.GetStatistics( True, True ) print "[ STATS ] = Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( stats[0], stats[1], stats[2], stats[3] ) [ STATS ] = Minimum=17.000, Maximum=255.000, Mean=220.586, StdDev=39.705 [ STATS ] = Minimum=64.000, Maximum=255.000, Mean=214.975, StdDev=36.926 [ STATS ] = Minimum=45.000, Maximum=255.000, Mean=179.029, StdDev=68.234

But if you use band.ReadAsArray() (= Numpy array)

for j in range(bands): band = f.GetRasterBand(j+1) data = band.ReadAsArray() print "[ Numpy ] = Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % (data.min(), data.max(), data.mean(), data.std()) [ Numpy ] = Minimum=0.000, Maximum=255.000, Mean=220.477, StdDev=42.584 [ Numpy ] = Minimum=31.000, Maximum=255.000, Mean=214.955, StdDev=39.558 [ Numpy ] = Minimum=0.000, Maximum=255.000, Mean=178.856, StdDev=69.535

Why? The problem is (GDAL does not ignore NoData value)

GetStatistics will reuse previously computed statistics if they exist (i.e computed before you set the NoData value). You can use stats = band.ComputeStatistics(0) instead of GetStatistics to force the statistics to be recomputed.

for j in range(bands): band = f.GetRasterBand(j+1) stats = band.ComputeStatistics(0) print "[ STATS ] = Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( stats[0], stats[1], stats[2], stats[3] ) [ STATS ] = Minimum=0.000, Maximum=255.000, Mean=220.477, StdDev=42.584 [ STATS ] = Minimum=31.000, Maximum=255.000, Mean=214.955, StdDev=39.558 [ STATS ] = Minimum=0.000, Maximum=255.000, Mean=178.856, StdDev=69.535
Improve this answer
2

...Or you could just convert it to an ESRI Ascii Raster and achieve effectively the same result in much less time.

Here's an example of an ascii raster from the documentation:

ncols 480 nrows 450 xllcorner 378923 yllcorner 4072345 cellsize 30 nodata_value -32768 43 2 45 7 3 56 2 5 23 65 34 6 32 54 57 34 2 2 54 6 35 45 65 34 2 6 78 4 2 6 89 3 2 7 45 23 5 8 4 1 62 ...

GDAL can do the conversion very fast, you can then read the file, or whatever is required. There is nothing wrong with the other answer. I only suggest this because I find NUMPY very slow for cell-by-cell operations.

Improve this answer
3
  • Could you elaborate? Commented Dec 8, 2016 at 17:28
  • This is another option, but I preferred to adapt the answer starting from its code =) Commented Dec 8, 2016 at 17:34
  • Yeah - your answer is perfectly fine/correct. I'm just used to processing way too much data and look for alternative routes that save time (though are often admittedly a little left-field) Commented Dec 8, 2016 at 17:36
0

I would like to create a new raster that contains the sum of all pixelvalues of the raster contained in a folder. Is it possible?

os.chdir(r'C:/TifFolder') li_rasters = [raster for raster in os.listdir(os.getcwd()) if os.path.splitext(raster)[-1] == '.tiff'] #print (li_rasters) #final_band=(r'C:/Users/KIFF/Desktop/These/data/Results/result.tiff') for raster in li_rasters: ......
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.