2

I wrote a function to determine the proportion of raster pixels that contain each of a list of values that I would be interested in. I then wrote a function to apply the previous function to a 3d numpy array. I would like to know how to improve this and also whether it would be easier to apply a function like this to a list of raster files, rather than create a numpy array.

A test array and sample values:

import numpy as np
test_arr = np.random.randint(0, 200, 200).reshape(2,10,10)
values = [test_arr[1][1][1], test_arr[0][0][0]]

Here are my functions:

def pixel_props(rast, values):
"""this function outputs a list with the proportions of 
the total raster pixels that contain each value."""
prop_list = []
size = float(rast.size)
for i in values:
    temp_count = np.count_nonzero(rast[rast == i])
    prop_list.append(temp_count/size)
return prop_list

def raster_props(array, no_data):
"""this will apply the pixel_props function to a bunch of layers in a numpy 3d array"""
    master_matrix = [[]]
    master_matrix.append(no_data)
    for layer in array:
        master_matrix.append(pixel_props(layer, no_data))
    return master_matrix
  • Do the functions produce the expected outcome - if so which part do you need to improve specifically? – Kersten Oct 19 '15 at 14:44
  • @Kersten the function produces a list of lists, which I believe is what I want. I wondered if it would be preferable to write a function that iterates over files in a file list without using a numpy array. I also wondered if there were ways to make the script more concise, i.e. do the same using only one function. – RyanM Oct 19 '15 at 16:25
  • @Kersten, I would also like to do this for a large number of large rasters, and do not have the memory to create the numpy 3d arrays. If you or others could propose an alternate solution, that would be helpful. – RyanM Oct 19 '15 at 18:36
1

If you are running into memory issues it is indeed a good idea to loop over the individual raster files instead of reading all at once.

Here is how you would apply your pixel_props function to a directory of raster images by reading one after the other with rasterio. This assumes each image is a single-layer GeoTiff in a folder.

import os
import glob
import rasterio

def pixel_props(rast, values):
    """this function outputs a list with the proportions of 
    the total raster pixels that contain each value."""
    prop_list = []
    size = float(rast.size)
    for value in values:
        count = np.sum(rast == value)
        prop_list.append(count/size)
    return prop_list

# create a list for the location of each GeoTiff
tif_dir = "/geo/tiff/folder"
tif_list = sorted(glob.glob(os.path.join(tif_dir, "*.tif")))

# list of values you are interested in
val_list = [1, 2, 3]

# loop over the images and calculate the pixel_props
results = []
for tif_file in tif_list:
    with rasterio.open(tif_file, 'r') as tif:
        tif_arr = tif.read()
        results.append(pixel_props(tif_arr, val_list))

Since you have also been asking in another thread how to plot these results, this is an example for a MODIS quality layer time-series, where 1 means marginal data, 2 means snow and ice and 3 is cloud cover.

from matplotlib import pyplot as plt
plt.figure(figsize=(11,7))
for i, val in enumerate(val_list):
    plt.plot(np.array(results)[:,i], label="value %s" % val)
plt.legend(loc="upper right")
plt.show()

enter image description here

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.