1

I'm working on the below script. What I'm trying to do is create a stats.txt file output which has the Land Cover Class ID, Fractional Coverage (of that class) and the number of pixels (of that class) in a three column .txt so that I can pass the file, along with another .tif file I've created into the next stage of my classification project.

I've muddled through a set of errors and have solved a few (see: size-1 array problem).

My current problem is two-fold, my output file just creates a string of numbers not separated by anything so I'm not sure what that is caused by, but I've tried various methods to write to a file (as you can see below).

My main issue is that when I run the code in part, in ipython, I come across the issue that my input .tif file (a CEH landcover map product: https://eip.ceh.ac.uk/lcm/lcmdata) when read, is producing an array of just zeros and I'm not sure that is right.

In [3]: print lc
[[[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 ...

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]]

Below is the main bulk of the code I have thus far.

How do I fix this problem?

#!/usr/bin/env python

import sys
import calendar
import os

from osgeo import gdal
import numpy as np
from scipy.stats import mode

from IPython import embed

GDAL2NUMPY = {  gdal.GDT_Byte      :   np.uint8,
                gdal.GDT_UInt16    :   np.uint16,
                gdal.GDT_Int16     :   np.int16,
                gdal.GDT_UInt32    :   np.uint32,
                gdal.GDT_Int32     :   np.int32,
                gdal.GDT_Float32   :   np.float32,
                gdal.GDT_Float64   :   np.float64,
                gdal.GDT_CInt16    :   np.complex64,
                gdal.GDT_CInt32    :   np.complex64,
                gdal.GDT_CFloat32  :   np.complex64,
                gdal.GDT_CFloat64  :   np.complex128
              }


#Open the training LC map file.

fname = sys.argv[1]
lc_dataset = gdal.Open(fname)
lc = lc_dataset.ReadAsArray()
#lc = np.array(lc)

#Calculating total number of pixels with a valid Land Cover ID.

fill_value = 0
number_of_pixels = np.where(lc != fill_value)[0].shape[0]

#Get the number of classes and corresponding IDs.

lc_classes = np.unique(lc)

#Split each class into its contituante pixel and write result to file.

f = open("LCM2015_stats_test.txt", "w+")

for classID in lc_classes:
    lc_class_pixels = np.where(lc == classID)[0].shape[0]
    FractionalCover = lc_class_pixels/number_of_pixels
    DAT = np.column_stack((classID, FractionalCover, lc_class_pixels))
    np.savetxt('LCM2015_stats_text_2.txt', DAT, delimiter = "", fmt="%s")
    #f.write(classID, FractionalCover, lc_class_pixels)
    break
f.close()

If I can provide any extra information please let me know.

3

Firstly, why do you think the dataset is being read as all 0's? If you're just looking at the output of print(lc) then you're not seeing the whole array just snippets from the edges (3 elements from each end of each dimension by default). In the example below, I print out the array and it just shows the zeroes at each edge of the array, but then I print the max class value and the class count and it shows that there's more data in there than shown by the print(lc).

Secondly, you're not supplying any delimiters when writing out the file so no delimiters get written. The example below uses the csv module, but np.savetxt has a delimiters argument or you could just join the values into a delimited string yourself f.write(','.join([classID, FractionalCover, lc_class_pixels])+'\n') etc...

import csv 
from osgeo import gdal
import numpy as np

ds = gdal.OpenEx('lcm2015_ni_1km_dominant_target_class.img')
lc = ds.ReadAsArray()

fill_value = 0
number_of_pixels = np.where(lc != fill_value)[0].shape[0]

#Get the number of classes and corresponding IDs.
lc_classes = np.unique(lc)

#Split each class into its constituent pixel and write result to file.
with open("LCM2015_stats_test.csv", "w") as f:

    writer = csv.writer(f)
    writer.writerow(['classID', 'FractionalCover', 'lc_class_pixels']) # Header

    for classID in lc_classes:
        lc_class_pixels = np.where(lc == classID)[0].shape[0]
        FractionalCover = lc_class_pixels/number_of_pixels
        DAT = np.column_stack((classID, FractionalCover, lc_class_pixels))
        #np.savetxt('LCM2015_stats_text_2.txt', DAT, delimiter = "", fmt="%s")
        writer.writerow([classID, FractionalCover, lc_class_pixels])

print('array:\n', lc, '\n')
print('max:', lc.max(), '\n')
print('number_of_pixels:\n', number_of_pixels, '\n')
print('lc_classes:', lc_classes, '\n')

!cat LCM2015_stats_test.csv  # The "!" is IPython syntax to run a shell command, doesn't work in standard python interpreter/scripts

array:
 [[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]] 

max: 21 

number_of_pixels:
 14865 

lc_classes: [ 0  1  2  3  4  5  6  7  9 10 11 12 13 14 15 16 17 18 20 21] 

classID,FractionalCover,lc_class_pixels
0,1.0450723175243861,15535
1,0.007198116380760175,107
2,0.04803229061553986,714
3,0.036057854019508916,536
4,0.6708375378405651,9972
5,0.007332660612176253,109
6,6.727211570803902e-05,1
7,0.02536158762193071,377
9,0.010023545240497814,149
10,0.03477968382105617,517
11,0.057315842583249246,852
12,6.727211570803902e-05,1
13,0.015472586612848975,230
14,0.041439623276152034,616
15,0.0003363605785401951,5
16,0.003363605785401951,50
17,0.0013454423141607804,20
18,0.006054490413723511,90
20,0.003901782711066263,58
21,0.031012445341405986,461

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.