Mailman 3 July 2010 - NumPy-Discussion

Datarray bug tracker
by Keith Goodman July 24, 2010

July 24, 2010

Report datarray bugs here: http://github.com/fperez/datarray/issues A datarray is a subclass of a Numpy array that adds the ability to label the axes and to label the elements along each axis.

1 0

reviewers needed for NumPy
by Joe Harrington July 23, 2010

July 23, 2010

Hi folks, We are (finally) about to begin reviewing and proofing the NumPy docstrings! This is the final step in producing professional-level docs for NumPy. What we need now are people willing to review docs. There are two types of reviewers: Technical reviewers should be developers or *very* experienced NumPy users. Technical review entails checking the source code (it's available on a click in the doc wiki) and reading the doc to ensure that the signature and description are both correct and complete. Presentation reviewers need to be modestly experienced with NumPy, and should have some experience either in technical writing or as educators. Their job is to make sure the docstring is understandable to the target audience (one level below the expected user of that item), including appropriate examples and references. Review entails reading each page, checking that it meets the review standards, and either approving it or saying how it doesn't meet them. All this takes place on the doc wiki, so the mechanics are easy. Please post a message on scipy-dev if you are interested in becoming a reviewer, or if you have questions about reviewing. As a volunteer reviewer, you can put as much or as little time into this as you like. Thanks! --jh-- for the SciPy Documentation Project team

1 0

Audiolab 0.11.0
by David Cournapeau July 23, 2010

July 23, 2010

Hi, I am pleased to announce the 0.11.0 release of audiolab scikits, the package to read/write audio file formats into numpy. This release has barely no changes compared to 0.10.x series, but it finally fixes annoying windows issues (which ended up being mingw bugs). Source tarball and python 2.6 windows installer are available on pypi, cheers, David

1 0

Control format for array of integers
by John Reid July 23, 2010

July 23, 2010

Hi, My array of integers is printed like this by default in numpy: array([[ 4.75000000e+02, 9.50000000e+02, -1.00000000e+00], [ 2.60900000e+03, 9.50000000e+02, -7.09000000e+02]]) Can I set an option so that numpy never uses this scientific notation or raise the threshold at which it kicks in? Thanks, John.

2 1

Finding Unique Pixel Values
by Ian Mallett July 23, 2010

July 23, 2010

Hi, So, I'm working on a radiosity renderer, and it's basically finished. I'm now trying to optimize it. Currently, by far the most computationally expensive operation is visibility testing, where pixels are counted by the type of patch that was drawn on them. Here's my current code, which I'm fairly sure can be optimized: line 1) pixel_data = np.cast['int16'](255.0*glReadPixels(0,0,patch_view.width,patch_view.height,GL_RGB,GL_FLOAT)).reshape((-1,3)) line 2) pixel_data[:,1] *= 255 line 3) pixel_data[:,0] *= (255*255) line 4) pixel_data = np.sum(pixel_data,1) line 5) for pixel in pixel_data: line 6) if pixel != 0: line 7) try: self.visibility[poly1num][pixel-1] line 8) except: self.visibility[poly1num][pixel-1] = 0 line 9) self.visibility[poly1num][pixel-1] += 1 The basic idea: Convert unique an array of RGB triplets into and array of unique numbers. The triplet [0,0,1] maps to 1 and the triplet [255,255,255] maps to 16581375. The triplet [0,0,0], mapping to 0, is always discarded. The resulting array of integers in the range [1,16581375] are the (index+1)'s of a list, self.visibility[poly1num]. For each integer "i" in the array, the value in self.visibility[poly1num][i-1] must be incremented by 1. By line: line 1: convert the read pixels from the range [0.0,1.0] (as floats) to the range [0,255] (as integers). Change the shape from (width,height,3) to (width*height,3). lines 2, 3, and 4: the pixels are an index to patches. The green channel counts 255 times as much. The red, 255 times as much as the green. This way, 255^3 possible pixel types can be had. Sum the red, green, and blue channels to arrive at a single number, which is the patch index in the range [0,255^3]. line 5: for each index in the data line 6: I'm reserving the index [0,0,0] for no patch, so if the index is that color, no patch was drawn. lines 7, 8, and 9: increment self.visibility[poly1num][pixel-1] by 1. I have a feeling that this code is likely the bottleneck. It's taking about an hour to execute this process 8720 times when width and height are both 512. I have a feeling my code can be greatly improved in many ways--if not all over, certainly in lines 5 through 9. Help? Thanks, Ian

3 5

Re: [Numpy-discussion] Crosstabulation
by Ionut Sandric July 22, 2010

July 22, 2010

Thank you Zack: By raster data I mean classified slope gradient (derived from a dem), landuse-landcover, lithology etc. A crosstabulation analysis will give me a table with the common areas for each class from each raster and this will go into other analysis. I can do it with other softwares (like ArcGIS DEsktop etc), but I would like to have all with numpy or to build something on top of numpy Thanks's again Ionut ----- Original Message ----- From: "Zachary Pincus" <zachary.pincus(a)yale.edu> To: "Discussion of Numerical Python" <numpy-discussion(a)scipy.org> Sent: Wednesday, July 14, 2010 9:42:49 PM GMT +02:00 Athens, Beirut, Bucharest, Istanbul Subject: Re: [Numpy-discussion] Crosstabulation Hi Ionut, Check out the "tabular" package: http://parsemydata.com/tabular/index.html It seems to be basically what you want... it does "pivot tables" (aka crosstabulation), it's built on top of numpy, and has simple data IO tools. Also check out this discussion on "pivot tables" from the numpy list a while ago: http://mail.scipy.org/pipermail/numpy-discussion/2007-August/028739.html One question -- what do you mean by "raster data"? In my arena, that usually means "images"... and I'm not sure what a crosstabulation on image data would mean! Zach On Jul 14, 2010, at 2:28 PM, Ionut Sandric wrote: > > Sorry, the first email was sent before to finish it... > > > Hi: > > I have two raster data and I would like to do a crosstabulation > between them and export the results to a table in a text file. Is it > possible to do it with NumPy? Does someone have an example? > > Thank you, > > Ionut > > > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion(a)scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion(a)scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

6 9

Custom ufuncs with Axis argument
by John Salvatier July 22, 2010

July 22, 2010

What is the easiest way to give a custom ufunc an axis argument? I have looked around the UFunc API, but I have not seen anything related to this.

1 1

list to array is slow
by marco cammarata July 22, 2010

July 22, 2010

Hi, any idea why the simple code below is so slow ? import numpy as n from time import time as t dims = (640,480) m = n.random.random( dims ) l=[] for i in range(200): l.append(m) t0=t() b=n.array(l) print t()-t0 To convert the list into an array takes about 5 sec ... Thanks, marco

3 2

summarizing blocks of an array using a moving window
by Robin Kraft July 22, 2010

July 22, 2010

Hello all, The short version: For a given NxN array, is there an efficient way to use a moving window to collect a summary statistic on a chunk of the array, and insert it into another array? The long version: I am trying to resample an image loaded with GDAL into an NxN array. Note that this is for statistical purposes, so image quality doesn't matter. For the curious, the image is derived from satellite imagery and displays a map of hotspots of tropical deforestation at 500m resolution. I need to assign a count of these deforestation 'hits' to each pixel in an existing array of 1000m pixels. Let's say the image looks like this: np.random.randint(0,2, 16).reshape(4,4) array([[0, 0, 0, 1], [0, 0, 1, 1], [1, 1, 0, 0], [0, 0, 0, 0]]) I want to use a square, non-overlapping moving window for resampling, so that I get a count of all the 1's in each 2x2 window. 0, 0, 0, 1 0, 0, 1, 1 0 3 => 2 0 1, 1, 0, 0 0, 0, 0, 0 In another situation with similar data I'll need the average, or the maximum value, etc.. My brute-force method is to loop through the rows and columns to get little chunks of data to process. But must be a more elegant way to do this - it's probably obvious too ... (inelegant way further down). Another way to do it would be to use np.tile to create an array called "arr" filed with blocks of [[0,1],[2,3]]. I could then use something like this to get the stats I want: d0 = img[np.where(arr==0)] d1 = img[np.where(arr==1)] d2 = img[np.where(arr==2)] d3 = img[np.where(arr==3)] img_out = d0 + d1 + d2 + d3 This would be quite snappy if I could create arr efficiently. Unfortunately np.tile does something akin to np.hstack to create this array, so it isn't square and can't be reshaped appropriately (np.tile(np.arange(2**2).reshape(2, 2), 4)): array([[0, 1, 0, 1, 0, 1, 0, 1], [2, 3, 2, 3, 2, 3, 2, 3]]) Inefficient sample code below. Advice greatly appreciated! -Robin import numpy as np from math import sqrt from time import time def rs(img_size=16, window_size=2): w = window_size # making sure the numbers work if img_size % sqrt(img_size) > 0: print "Square root of image size must be an integer." print "Sqrt =", sqrt(img_size) return elif (img_size / sqrt(img_size)) % w > 0: print "Square root of image size must be evenly divisible by", w print "Sqrt =", sqrt(img_size) print sqrt(img_size), "/", w, "=", sqrt(img_size)/w return else: rows = int(sqrt(img_size)) cols = int(sqrt(img_size)) # create fake data: ones and zeros img = np.random.randint(0, 2, img_size).reshape(rows, cols) # create empty array for resampled data img_out = np.zeros((rows/w, cols/w), dtype=np.int) t = time() # retreive blocks of pixels in img # insert summary into img_out for r in xrange(0, rows, w): # skip rows based on window size for c in xrange(0, cols, w): # skip columns based on window size # calculate the sum of the values in moving window, #insert value into corresponding pixel in img_out img_out[r/w, c/w] = np.int(np.sum(img[r:r+w, c:c+w])) t1= time() elapsed = t1-t print "img shape:", img.shape print img, "\n" print "img_out shape:", img_out.shape print img_out print "%.7f seconds" % elapsed rs(imgage_size=16, window=2) rs(81, 3) rs(1000000) #rs(4800*4800) # very slow

5 6

Re: [Numpy-discussion] summarizing blocks of an array using a moving window
by Robin Kraft July 22, 2010

July 22, 2010

Vincent, Pauli, > From: Vincent Schut <schut(a)sarvision.nl> > - an other option would be some smart reshaping, which finally gives you > a [y//2, x//2, 2, 2] array, which you could then reduce to calculate > stats (mean, std, etc) on the last two axes. I *think* you'd have to > first reshape both x and y axes, and then reposition the axes. An > example: > a = gdal_array.BandReadAsArray(blabla) > y,x = a.shape #y and x need be divideable by 2! > b = a.reshape(y/2, 2, x/2, x).transpose(0,2,1,3).reshape(y/2, x/2, 4) > bMean = b.mean(axis=-1) > bMax = ......etc. You and Pauli agree on this - seems like a good option. > - a third option would be to create an index array, which has a unique > value per 2x2 square, and then use histogram2d. This would, if you use > its 'weight' functionality, at least enable you to get efficient counts > and sums/means. Other stats might be hard, though. Hmmmm I don't get this, but I'll experiment. I've seen some really useful things done with histogram2d so it seems worthwhile to figure out. > Message: 6 > From: Pauli Virtanen <pav(a)iki.fi> > > Let's say the image looks like this: np.random.randint(0,2, > > 16).reshape(4,4) > > > > array([[0, 0, 0, 1], > > [0, 0, 1, 1], > > [1, 1, 0, 0], > > [0, 0, 0, 0]]) > > > > I want to use a square, non-overlapping moving window for resampling, so > > that I get a count of all the 1's in each 2x2 window. > > > > 0, 0, 0, 1 > > 0, 0, 1, 1 0 3 > > => 2 0 > > 1, 1, 0, 0 > > 0, 0, 0, 0 > > > > In another situation with similar data I'll need the average, or the > > maximum value, etc.. > > Non-overlapping windows can be done by reshaping: > > x = np.array([[0, 0, 0, 1, 1, 1], > [0, 0, 1, 1, 0, 0], > [1, 1, 0, 0, 1, 1], > [0, 0, 0, 0, 1, 1], > [1, 0, 1, 0, 1, 1], > [0, 0, 1, 0, 0, 0]]) > > y = x.reshape(3,2,3,2) > y2 = y.sum(axis=3).sum(axis=1) This is perfect - and so fast! Thanks! Now I just have to understand why it works ... Can anyone recommend a tutorial on working with (slicing, reshaping, etc.) multi-dimensional arrays? Stefan's slides are beautiful but my brain starts to hurt if I try to follow them line by line. -Robin

1 0