Mailman 3 April 2008 - SciPy-User

Re: [SciPy-user] [AstroPy] FITS images with header-supplied axes?
by James Turner April 2, 2008

April 2, 2008

Hi everyone, I agree that it would be good for the community to get into the habit of exposing work on astronomical Python modules via SciPy (for example). It seems that the obvious place for this is "AstroLib", which has been part of the SciPy site for some time, but doesn't seem to have gathered a lot of momentum yet: http://scipy.org/AstroLib I'm also interested in WCS, as critical functionality that hasn't quite crystallized yet. Wrapping WCSLib did seem like an obvious direction when I first looked into this and talked to Perry, but it seemed that there may be licensing issues to investigate, with it being GPL (which I personally like, but I'm not sure about internal policy on relicensing everything). At Gemini, one thing we have been working on is a data access class that takes care of some logistics like figuring out the "type" (meaning instrument, for example, not numerical type) of data, accessing standard header information in an instrument-agnostic way, iterating over things like nod sets, FITS extensions and so on. In future, the class may also support array operations with automatic propagation of things like variance and header info. This is still really alpha code at the moment though and I'm not announcing it officially. Later in the year, we will be working on some generic spectroscopic data reduction tasks. These things are planned to operate both in a pipeline context and interactively. I also have a little script for mosaicing IFU data that might be of interest to someone, but it still needs converting to NumPy. That will probably work its way into the official Gemini DR package at some point. By the way, those of us at AURA sites will be meeting in June to discuss collaboration on new DR infrastructure (in general). Cheers, James. -- James E.H. Turner Gemini Observatory Southern Operations Centre, Casilla 603, Tel. (+56) 51 205609 La Serena, Chile. Fax. (+56) 51 205650

2 1

Re: [SciPy-user] FITS images with header-supplied axes?
by Perry Greenfield April 2, 2008

April 2, 2008

On Mar 28, 2008, at 8:28 PM, Keflavich wrote: > Is there any plotting routine in scipy / matplotlib that can plot a > fits image with correct WCS coordinates on the axes? I know pyfits > can load fits files, astLib has routines to interpret header coordinates, and I think you can make the axes different using > matplotlib transforms, but is there anything that puts all three > together currently available? > Thanks, > Adam Well, we (STScI) recently wrapped WCSLIB to obtain a mapping function between pixel and sky coordinates for python (you can find it as pywcs in astrolib on scipy; that may have been what you were referring to). But I'm not sure you understand what you are asking for with regard to matplotlib. The new transforms stuff should make it much easier to display the sky coordinates in the interactive display. The axis labeling is a different matter. Suppose your image (let's say it's 1Kx1K for the sake of discussion) is rotated 45 degrees with regard to north (either way, it doesn't really matter). What would you expect to see for axis labels? I don't think it is at all obvious how people would want labeling to be done along the edges of the image. I can imagine someone wanting axes or grids superimposed on the image itself, but that's not quite the same thing. Do you want the image rotated so that it is resampled on to RA and Dec and displayed that way? In any event, no we haven't yet done anything to try to integrate all three things. Among other things we wanted to make sure that the api for the wcs info was suitable before doing a lot with it (and in the meantime, Mike is working on rewriting drizzle which is taking a lot of his time). Perry

2 3

kmeans2 bug for 1D data?
by Lin Shao April 2, 2008

April 2, 2008

I understand that kmeans2 is usually used for multidimensional vector space, but it is sometimes useful for 1D clustering, such as clustering the pixels of an image based solely on pixel intensities. And kmeans2 does theoretically support 1D, as stated in its function doc: def kmeans2(data, k, iter = 10, thresh = 1e-5, minit = 'random', missing = 'warn'): """ ...... data : ndarray Expect a rank 1 or 2 array. Rank 1 are assumed to describe one dimensional data, rank 2 multidimensional data, in which case one row is one observation. ........ But the truth is if data is 1D and if minit is 'random', there's an error message when calling vq(data, code) in _kmeans(), because apparently code is a rank-2 array. The cause is in _krandinit() where the return value, x, is rank-2 no matter input data is rank 1 or 2. I "fixed" it by replacing "return x" with "return x.squeeze()". I'm not sure if that's the right way. -- The magic of the microscope is not that it makes little creatures larger, but that it makes a large one smaller.

1 0

Astronomical Python Projects (was Re: [AstroPy] FITS images with header-supplied axes?)
by Tom Bridgman April 2, 2008

April 2, 2008

Since I'm at home now, I'll describe my 'recreational' python projects. I'm probably the really weird one on the list. ;^) My major use is in developing simple models/simulations & data analysis exercises which can be used in debunking crank astronomy & physics (see my .sig below). I'm using scipy/numpy/pyfits/PIL/PyX for generating graphics from simple models. Here's some recent examples: - electron-proton-neutron gas equation of state (neutron stars can't exist claims) - power spectra of simple cosmological models (redshift quantization claims) - simple galactic supernova distribution models (not enough SNR claims) - simple radiation transport simulations - simple E&M field configuration simulations (trying out FiPy for this) - simulated 'observations' based on various real & crackpot cosmological models - nuclear decay & nuclear reaction networks. For this project, I don't have a strong need for precision tools. I don't use PyRAF, or the more complex tools available. I want most of my projects to be easily reproducible by students and instructors at various levels. FV & DS9 have sufficed for viewing images generated by image-based models to unprocessed Hubble imagery. The frameworks for some of these might be of use to others, but I'm not sure the best way to package them. I've looked at some of the data-table readers promoted on the list, but have generally found them overkill for my needs (hence I re-invent the wheel for much of this). More recently, I've been working with the pyPOV wrapper to generate nice visualizations from models using the POVray ray-tracer. I do have an interest in a good implementation of WCS in case I want to overlay a model with actual data. I have a simple UTC time class but it is *not* a precision object (unsuitable for say, pulsar or x- ray timing observations). I'll post the work list tomorrow. Thanks, Tom On Apr 1, 2008, at 4:32 PM, Russell E Owen wrote: > At 1:51 PM -0400 2008-04-01, Bridgman, William T. wrote: >> >> Would there be any interest in members of the list publishing a short >> description of what types of modules they are designing in their own >> work? It might be worthwhile for coordination & possible >> collaborations. My requirements for work projects are quite >> different from my recreational & educational python projects. > > Good idea. > > I'm part of the team working on data processing for the LSST. This > number crunching is done in C++ and the high level operations are > done in python. This work will probably be primarily of interest to > those processing a lot of data because it is pipeline-oriented. > > The final product will include: > - basic types for images, masks and masked images > - image linearization (bias subtraction, flat fielding, etc.) > - wcs determination > - source detection > - image subtraction > > We intend it to be usable on a variety of data sources (not just > LSST) since that's the only way to test it. Nonetheless it will > probably take a bit of work to massage the data headers. > > As far as wcs goes: right now we use wcslib and have a limited python > interface on it. I suspect if anything better came along we'd be > happy to switch. (at least at the C++ level; I'm not sure about the > Python level). > > -- Russell > _______________________________________________ > AstroPy mailing list > AstroPy(a)scipy.org > http://lists.astropy.scipy.org/mailman/listinfo/astropy -- Dealing with Creationism in Astronomy http://homepage.mac.com/cygnusx1 cygnusx1(a)mac.com "They're trained to believe, not to know. Belief can be manipulated. Only knowledge is dangerous." --Frank Herbert, "Dune Messiah"

1 0

paths relative to a python package
by Zane Selvans April 2, 2008

April 2, 2008

I know this isn't scipy specific but... If you have a data file that's distributed with a python package, how do you refer to it relative to the location that the package ultimately gets installed? Is there a package_root variable or something? e.g. in my package directory let's say I have: __init__.py MyModule.py datadir and within datadir: data1.dat data2.dat ... How do I, within MyModule.py tell the program to open one of the data files? If I just do open('datadir/data1.dat') it'll try and open that dir/file relative to wherever the python program that imported MyModule is running, right? Which isn't what I want. Thanks for any insight, Zane -- Zane Selvans Amateur Human zane(a)ideotrope.org 303/815-6866 PGP Key: 55E0815F

2 2

scipy and blitz
by nicky van foreest April 1, 2008

April 1, 2008

Dear all, After running laplace.py (a part of the scipy performance tutorial) I get the following result: chuck~/tmp/perfpy> python laplace.py Doing 100 iterations on a 500x500 grid numeric took 3.35 seconds blitz took 8.61 seconds inline took 7.84 seconds fastinline took 0.75 seconds fortran took 0.42 seconds pyrex took 0.73 seconds whereas the results on http://www.scipy.org/PerformancePython give this: Numeric 29.3 Blitz 9.5 Inline 4.3 Fast Inline 2.3 Python/Fortran 2.9 Pyrex 2.5 Does anybody have clue about why the blitz code works slower on my machine than numpy? I searched the web on what might be a reason and came across a mail by Travis Oliphant. He also experiences relatively slow behavior of blitz, but after recompiling python with debugging info blitz worked well again. Thus, I checked my installation of python, but it is compiled without debugging info. As a further check I contacted the scipy maintainer of gentoo (see cc). He gets the same relative performance as I do. Thanks for any info. bye Nicky

2 4

kmeans2 random initialization
by Alec Solway April 1, 2008

April 1, 2008

Hello, In using kmeans2 with the "random" initialization method, I often get: numpy.linalg.linalg.LinAlgError: Matrix is not positive definite - Cholesky decomposition cannot be computed I haven't dug too deeply into the implementation; I'm wondering if there's a gotcha someone can point out quickly here. What exactly is the "random" initialization method doing with the data set, and what constraints must the data have so as to keep whatever intermediate matrix it computes to be positive definite? Thanks.

3 3

Impulse/Step response troubles
by Kumar Appaiah April 1, 2008

April 1, 2008

Dear SciPy users, (Long mail) I am having several problems in working out the impulse response of a simple function in SciPy. I would request you to please point out whether the fault lies with SciPy or me. The question is to determine (plot) the step and ramp responses of K / (s^2 + 8s + K) for K = 7, 16 and 80 So, the following code is used, and I'll keep changing a line: <code> from scipy import * from pylab import * K = 16.0 # should be redone with 7.0, 16.0 and 80.0 r = signal.impulse(([K],[1.0,8.0,K]), T=r_[0:5:0.001]) plot(r[0], r[1], linewidth=2) show() </code> The above code plots the impulse response of the given function. Now, on to the action: 1. Running the above code with K = 7.0 and K = 16.0 gives expected results. However, running the code with K = 16.0 doesn't work; if K = 16.0, the response should be 16 * t * exp(-4t) * u(t), which it isn't. Changing K to 16 + or - .00000000001 fixes it. There surely is some problem when the value of K is such that a double root is hit. 2. For the step response, I change "impulse" in the above code to "step", things seem all right. Again, for K = 16, it doesn't work. I then switch back to impulse, but make the second (denominator) array [1.0, 8.0, K, 0], and still the same thing is observed. 3. To get the ramp response, I now have the following code in line 8 of the above code. r = signal.step(([K],[1.0,8.0,K, 0]), T=r_[0:5:0.001]) But this causes the error: Traceback (most recent call last): File "<stdin>", line 8, in ? File "/usr/lib/python2.4/site-packages/scipy/signal/ltisys.py", line 498, in step vals = lsim(sys, U, T, X0=X0) File "/usr/lib/python2.4/site-packages/scipy/signal/ltisys.py", line 403, in lsim ATm1 = linalg.inv(AT) File "/usr/lib/python2.4/site-packages/scipy/linalg/basic.py", line 306, in inv if info>0: raise LinAlgError, "singular matrix" numpy.linalg.linalg.LinAlgError: singular matrix So let's switch back to impulse: r = signal.impulse(([K],[1.0,8.0,K, 0.0, 0.0]), T=r_[0:5:0.001]) Now, it seems to work, but I am yet to verify the correctness of the plots. Naturally, for K = 16.0, it doesn't work; I have to use a value close to it. I would really appreciate help in getting this debugged. Thanks a lot! Kumar -- Kumar Appaiah, 458, Jamuna Hostel, Indian Institute of Technology Madras, Chennai - 600 036

2 3

Re: [SciPy-user] FITS images with header-supplied axes?
by humufr＠yahoo.fr April 1, 2008

April 1, 2008

Hello, I just want to discuss about a problem that our python astronomer have and the answer from Perry is a very good example of it. There are too many projects which are doing exactly the same thing and essentially because nothing is centralized or worst, in the case of pywcs, not advertised. For example Perry told us about pywcs developped by STSCI (it was the first time I saw any reference to this project) but Adam spoke about astlib ( http://astlib.sourceforge.net/ ) which are different package with exactly the same goal and I did myself something similar (even if it was fast and dirtier). I think it's time to try to identified the most important task that the astronomer are needing and try to centralized all the effort at the same place. The astropy mail list is probably a good start as the astropy.scipy.org website. Perhaps we can start to identified the need and desiderata. We will have too many but that some most important or urgent must be identified. As example, we clearly needed the pywcs and thanks to STSCI we have it now. We also need a package to plot our data, images etc. Matplotlib is very good but does have a major problem (at least for my point of view), it's slow, very slow for big array or can't even produce the image if the image is too big. ( If I remember a precedent discussion, the problem is mainly due to Agg.) After we need to know, if possible, for which project if someone is doing something, who and how to contact him/her. So interested people can eventually help him/her to do it. It's seems that STSCI is doing most of the work but I'm pretty sure that other people are willingly to help them to extend python to be the ideal tool for astronomer. Just my 2 cents, Nicolas Le Sunday 30 March 2008 10:54:04 Perry Greenfield, vous avez écrit : > On Mar 28, 2008, at 8:28 PM, Keflavich wrote: > > Is there any plotting routine in scipy / matplotlib that can plot a > > fits image with correct WCS coordinates on the axes? I know pyfits > > can load fits files, astLib has routines to interpret header > > coordinates, and I think you can make the axes different using > > > matplotlib transforms, but is there anything that puts all three > > together currently available? > > > > Thanks, > > Adam > > Well, we (STScI) recently wrapped WCSLIB to obtain a mapping function > between pixel and sky coordinates for python (you can find it as pywcs in > astrolib on scipy; that may have been what you were referring to). > > But I'm not sure you understand what you are asking for with regard to > matplotlib. The new transforms stuff should make it much easier to display > the sky coordinates in the interactive display. The axis labeling is a > different matter. Suppose your image (let's say it's 1Kx1K for the sake of > discussion) is rotated 45 degrees with regard to north (either way, it > doesn't really matter). What would you expect to see for axis labels? I > don't think it is at all obvious how people would want labeling to be done > along the edges of the image. I can imagine someone wanting axes or grids > superimposed on the image itself, but that's not quite the same thing. Do > you want the image rotated so that it is resampled on to RA and Dec and > displayed that way? > > In any event, no we haven't yet done anything to try to integrate all three > things. Among other things we wanted to make sure that the api for the wcs > info was suitable before doing a lot with it (and in the meantime, Mike is > working on rewriting drizzle which is taking a lot of his time). > > Perry > _______________________________________________ > SciPy-user mailing list > SciPy-user(a)scipy.org > http://projects.scipy.org/mailman/listinfo/scipy-user

1 0