generating random selection of pixels in a raster then store pixel locations/values

[MargaretWooten]

Hi,

I have a large dataset of integer rasters, and was hoping to find a way to generate a "random" selection of pixels for four of the rasters, while making sure all 5 classes, which I've specified based off the pixel values, are equally represented with the selection. Each raster's values ranges from 0-200, with 0-100 representing percent tree coverage and 200 representing water (aka data I want to ignore for this case). I already added a new field to the raster called "TreeClass," and filled that field with numbers 1-6, where 1-5 are different levels/classes of tree coverage and 6 is water/no data values that I want to ignore. In all for each raster, I would like 10 pixels, 2 in each class 1-5.

In addition to randomly selecting these pixels, I would like to record their locations (plus values and classes, if possible), so that I can do a time series analysis on the selected pixels using the other rasters

Any suggestions on the best method for this? I have thought about several approaches:

1. Converting the four rasters to polygon based on TreeClass and then use the Create Random Points tool with the polygon as a constraint. However, after using the Raster to Polygon tool, it doesn't appear to have come out right, and is very slow. Additionally, I don't see a way in which this tool would allow me to store the pixel location (not the coordinates, but samples, lines location data of the pixel).

2. Converting the raster to numpy array and then using python tools to generate random pixels (not yet sure how I would go about this). That would make it easy to store the x, y location of the pixel and then use them to get the same pixel location's value for the other numpy arrays (arrays that I would convert from raster). However, I don't see a way to convert the raster to a structured array so that not only the pixel value would be present and accessible but the class value too. Perhaps I could convert to numpy array that stores the pixel value, somehow add a dimension to the numpy array, and then calculate the class value in python to store into the new dimension.

I'm hoping there's a simple solution and I'm just over-complicating things, as I tend to do. But I wanted to get some feedback that will help lead me in the right direction.

Thanks so much in advance.

[ChristopherBlinn1]

Margaret,

Are you familiar with Hawth's Analysis Tools?

There is a random selection tool that could help you.  You can read on it here.

To get points, you could convert your rasters to a point feature class using the Raster to Point tool.

When you create the points you could add two fields (type double, called lat/lng or X/Y, whichever is more appropriate), and calculate the geometry of each point.  This will add the location data to a field.  This will give you location and raster value for each point (essentially each pixel in your raster).

From there you could split the layer by attribute value (either manually or with a model), which will create separate point feature classes for each value in your TreeClass field.

Then you can run the random selection tool on each point feature class for each test.

Hope this helps!

Best,
Chris B.

[MargaretWooten]

Hi Chris,

Thanks for the suggestion! I looked at the Hawth's Analysis dialog and read that it is discontinued/shouldn't be used for later versions of Arc. I was thinking I use your instruction and convert the raster to point shapefile, add the lat/long fields, split the layers based on TreeClass then just use the 'Create Random Point' built in tool in Arc on each tree class shapefile.

However, working with these files (which are 4800x4800 pixels) as vectors is proving to be really slow. I successfully ran the Raster to Point on one raster, although it took forever, and splitting these up into 5 shapefiles based on class is taking even longer.

Do you have any other ideas that don't require converting to vector? Maybe something with python?

Thanks again!

[ChristopherBlinn1]

Margaret,

Off the top of my head, creating 10 random points then using the Extract Values to Points would work up until the point where 2 points must be in each class 1 thru 5.  The only other way I see to incorporate this approach is to create a loop that will create one random point at a time.  It will then check the TreeClass value after the Extract Values to Points function.  Then a counter variable (one for each class) will monitor which class the random point has.  Once the counter reaches 2, the point is not kept, and the loop continues to run and create 1 random point until each counter equals 2.  Once 10 points, 2 for each class are created, the files are merged into a feature class.  There's no telling how long this will take though.

I was able to get most features of Hawth's tools to work in ArcGIS 10 by getting the VBA extension installed as well.  If you're using 10.1, you could use GME.  This has incorporated many of the Hawth's functions, but requires R and does not run within ArcGIS.  I believe r.sample would be the best command to run.  This of course still relies on each TreeClass being it's own point feature class.

Sorry I can't be much more of a help.  The loop idea would work, but I'm not a python expert so I wouldn't be able to get you any solid sample code to start with.

Best,
Chris B.

[MargaretWooten]

Thanks again for your suggestion! I was able to successfully come up with a method for what I needed using RasterToNumPyArray and then doing everything else in Python. Here is my code, just for reference, and although I'm sure there is a more efficient way to go about this, it at least gave me what I needed:

import arcpy
import numpy
import random

VCtif = "D:\\Maggie Data\\VCF Time Series Data\\Trial\\phase2\\TC_2010.h21v02.tif"
VC = arcpy.RasterToNumPyArray(VCtif) # VC tif will be 1 of 4 VC tiffs, from loop

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# part a) generate random pixels
# create tuples for tree classes
c1 = numpy.where((VC >= 0) & (VC <= 10)) # class 1
c2 = numpy.where((VC >= 11) & (VC <= 20))
c3 = numpy.where((VC >= 21) & (VC <= 40))
c4 = numpy.where((VC >= 41) & (VC <= 60))
c5 = numpy.where((VC >= 61) & (VC <= 100))
classlist = [c1, c2, c3, c4, c5]

points = [] # empty list to store random points (will be a list of tuples)
# loop through c1...c5 and generate 2 random points for each
for c in classlist:
    X = c[0].tolist() # x list 
    Y = c[1].tolist() # y list

    x1= random.choice(X) # choose random point from class n
    ind = X.index(x1) # find location of this point
    y1 = Y[ind] # find corresponding y point
    points.append((x1, y1)) # add to points list

    x2= random.choice(X) # do the same for a second point in class n
    while x2 == x1:
        x2 = random.choice(X) # keep generating random second point to ensure it is different from x1
    ind2 = X.index(x2)
    y2 = Y[ind2]
    points.append((x2, y2))

##for i in range(0,10):
##    print VC[points] # ensure that values represent 5 class equally