tag:blogger.com,1999:blog-88012649367209393122024-03-13T10:27:42.713-06:00Python-AstroThis is a support for a lecture on Python given at the Instituto de Astronomia at the UNAM (Universidad Nacional Autonoma de Mexico) by Christophe Morisset.Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.comBlogger22125tag:blogger.com,1999:blog-8801264936720939312.post-57413657816412466612015-09-16T15:01:00.000-05:002015-09-16T15:01:10.649-05:00Links to cool pages on plotting with pythonA very interesting blog on python there: <a href="http://zulko.github.io/blog/2014/09/20/vector-animations-with-python/">http://zulko.github.io/blog/2014/09/20/vector-animations-with-python/</a><br />
<br />
A package to plot interactively: <a href="http://bokeh.pydata.org/en/latest/index.html">http://bokeh.pydata.org/en/latest/index.html</a><br />
<br />
And don't forget the impressive Mayavi package: <a href="http://mayavi.sourceforge.net/">http://mayavi.sourceforge.net/</a><br />
<br />
The Seaborn library is useful to make some nice plots: <a href="http://stanford.edu/~mwaskom/software/seaborn/">http://stanford.edu/~mwaskom/software/seaborn/ </a><br />
<br />
<br />
A kind of summary of some of these different solutions is made there: <a href="http://pbpython.com/visualization-tools-1.html">http://pbpython.com/visualization-tools-1.html</a><br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com19tag:blogger.com,1999:blog-8801264936720939312.post-87016737125189197312014-12-01T16:05:00.001-06:002014-12-01T16:05:23.189-06:00Nice book on python, numpy, ipython etc...A very up-to-date (august 2014) electronic book on python named "Introduction to Python for Econometrics, Statistics and Data Analysis" from Kevin Sheppard is available here:<br />
<a href="https://www.kevinsheppard.com/images/0/09/Python_introduction.pdf">https://www.kevinsheppard.com/images/0/09/Python_introduction.pdf</a><br />
It's mainly on econometric, but most of the tools described there are also useful for astronomers.<br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com1tag:blogger.com,1999:blog-8801264936720939312.post-25801958819291395042014-11-11T15:18:00.003-06:002016-05-31T15:35:53.481-05:00Sending requests to MySQL and receiving the result from python, using PyMySQL Modern astrophysics is using every day more big databases. One of the mostly used interface to databases is MySQl (or its recent free fork MariaDB). I present in this lecture a python library to deal with MySQL databases: PyMySQL. In the lecture the examples are using access to 3MdB <a href="https://sites.google.com/site/mexicanmillionmodels/">(https://sites.google.com/site/mexicanmillionmodels/</a>), which requires a password. You can ask for it to me, or adapt the example to connect to other databases. <br />
The lecture is here:<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_PyMySQL.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_PyMySQL.ipynb</a><br />
<br />
<br />
An introduction to MySQL can be found here: <a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/MySQL.pdf">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/MySQL.pdf</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-41212208285891476442014-10-08T21:31:00.003-05:002016-05-31T16:18:16.904-05:00Optimization, calling Fortran<h2>
2014 Python Lecture. Part IX</h2>
<br />
In this latest lecture of this series, I'll present some tools to optimize your code by CPU and memory profiling. It also contains some tips on using the python debugger.<br />
The notebook is there:<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Optimization.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Optimization.ipynb</a> <br />
<br />
<br />
I also give some indications on how one can call Fortran routines from within python, to accelerate the execution of some part the the code.<br />
Here are small examples: <a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Calling%20Fortran.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Calling%20Fortran.ipynb</a><br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-33652883910612213502014-10-01T13:41:00.000-05:002016-05-31T15:37:42.435-05:00Object Oriented Programing. Objects, classes, etc...<h2>
2014 Python Lecture. Part VIII</h2>
In this lecture I'll introduce the basic (and some not that basic) concepts of Object Oriented Programing. I'll use an example to show how to:<br />
<ul>
<li>use functions to do simple jobs</li>
<li>but use objects when things start to be more complex</li>
<li>define classes, objects, attributes, methods, etc...</li>
<li>use *args and **kwargs in functions calls</li>
<li>use the class variables</li>
<li>add functionalities to classes and objects</li>
<li>use class inheritance</li>
<li>use attributes properties</li>
</ul>
The notebook is here:<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/OOP.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/OOP.ipynb</a> Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-6856227164598164992014-09-25T13:47:00.000-05:002016-05-31T15:38:06.828-05:00The astropy library<h2>
2014 Python Lecture. Part VII</h2>
The Astropy Project is a community effort to develop a single core
package for Astronomy in Python and foster interoperability between
Python astronomy packages. More informations here: <a href="http://www.astropy.org/">http://www.astropy.org/</a><br />
<br />
In this lecture we will see some of the facilities of the astorpy library, including:<br />
<br />
<ul>
<li>Constants and Units</li>
<li>Data Table (a very useful one!)</li>
<li>Time and Dates</li>
<li>Etc...</li>
</ul>
<br />
The lecture is here:<br />
<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_astropy.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_astropy.ipynb</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-63948816326093221232014-09-18T13:03:00.001-05:002016-05-31T15:38:30.924-05:00Useful libraries<h2>
2014 Python Lecture. Part VI</h2>
<div>
This lecture will give some insights to the most useful python
libraries. It is NOT exhaustive, you have to read the corresponding
manual pages to find the best use you can have of them. The list of all
python-included libraries is here: <a href="https://docs.python.org/2/library/" target="_blank">https://docs.python.org/2/library/</a></div>
<div>
<br />
I mention in this lecture:<br />
<br />
<ul>
<li>time and datetime</li>
<li>timeit</li>
<li>os</li>
<li>sys</li>
<li>subprocess</li>
<li>glob</li>
<li>re</li>
<li>urllib2</li>
</ul>
</div>
<div>
The lecture is here:<br />
<a href="http://nbviewer.ipython.org/url/132.248.1.102/~morisset/Python_lectures/Useful_libraries.ipynb"> </a></div>
<div>
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Useful_libraries.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Useful_libraries.ipynb</a></div>
Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com1tag:blogger.com,1999:blog-8801264936720939312.post-72466814288752813042014-09-17T23:14:00.001-05:002016-05-31T15:39:18.196-05:00Introduction to Scipy<h2>
2014 Python Lecture. Part V</h2>
SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem
of open-source software for mathematics, science, and engineering. In
particular, these are some of the core packages: Numpy, Scipy library, Matplotlib, ipython, Simpy, and Panda.<br />
<br />
<br />
In this lecture, I will show exemples covering:<br />
<ul>
<li>Some useful methods</li>
<ul>
<li>nanmean</li>
<li>constants</li>
</ul>
<li>Integrations</li>
<li>Interpolations</li>
<li>2D-interpolations</li>
<li>data fitting</li>
<li>multivariate estimation</li>
</ul>
The lecture is here:<br />
<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Scipy.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Scipy.ipynb</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-90425094653026829562014-09-10T17:14:00.000-05:002016-05-31T15:39:44.277-05:00How to make plots, images, 3D, etc, using Matplotlib<h2>
2014 Python Lecture. Part IV</h2>
<br />
This lecture is dedicated to the plotting library matplotlib. The topics are:<br />
<br />
<br />
<ul>
<li>Simple plot</li>
<li>Controlling colors ans symbols</li>
<li>Overplot</li>
<li>Fixing axes limits</li>
<li>Labels, titles</li>
<li>Legends</li>
<li>The object oriented way to use Matplotlib</li>
<li>Scatter</li>
<li>log plots </li>
<li>Multiple plots</li>
<li>Everything is object</li>
<li>Error bars</li>
<li>Sharing axes</li>
<li>Histograms</li>
<li>Boxplots</li>
<li>Ticks, axes and spines</li>
<li>A plot inside a plot</li>
<li>Play with all the objects of a plot</li>
<li>Filled regions</li>
<li>2D-histograms</li>
<li>2D data sets and images</li>
<li>Contour</li>
<li>3D scatter plots</li>
<li>Saving plots</li>
<li>Access and clear the current figure and axe </li>
<li>What's happen when not in a Notebook? plt.show() and plt.ion() commands </li>
</ul>
The lecture Notebook is there:<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Matplotlib.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Matplotlib.ipynb</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-64107251489055823552014-09-03T16:02:00.003-05:002016-05-31T15:40:05.546-05:00Interacting with files: reading writing, ascii and fits<h2>
2014 Python lecture. Part III </h2>
<br />
It's time to play with files containing data! In this lecture, we'll see how to read and write files (ascii and fits).<br />
<br />
<ul>
<li>Reading a simple ASCII file</li>
<li>How to treat special rows (comments, header)</li>
<ul>
<li>classical way </li>
<li>using numpy.loadtxt</li>
<li>using numpy.genfromtxt</li>
</ul>
</ul>
<ul>
<li>Dealing with missing data</li>
<li>Data in a fixed size format</li>
<li>Writing files</li>
<ul>
<li>simple method</li>
</ul>
<li>Pickle files (python format)</li>
<li>FITS files</li>
</ul>
The ipython notebook is there:<br />
<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Interact%20with%20files.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Interact%20with%20files.ipynb</a><br />
<br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com8tag:blogger.com,1999:blog-8801264936720939312.post-11145147873561653472014-08-20T15:29:00.003-05:002016-05-31T15:41:15.015-05:00Introduction to Numpy<h2>
2014 Python lecture. Part II </h2>
<br />
The introduction to Numpy can be seen here:<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_numpy.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_numpy.ipynb</a><br />
<br />
The topics that are presented are:<br />
<br />
<ul>
<li>The Array class</li>
<ul>
<li>create an array</li>
<li>1D, 2D 3D arrays</li>
<li>creating array from scratch</li>
<li>arrays share memory (views)</li>
</ul>
<li>random generator</li>
<li>timing a command</li>
<li>slicing arrays</li>
<li>assignments</li>
<li>using masks</li>
<li>the where function</li>
<li>some operations with arrays</li>
<li>broadcasting</li>
<li>calling scripts</li>
<li>structured arrays and record arrays</li>
<li>NaN other ANSI values.</li>
</ul>
<div>
Any comments are welcome.</div>
<div>
Chris.Morisset a t Gmail.com</div>
<div>
<br /></div>
Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-18729640671057822942014-08-13T12:51:00.001-05:002016-05-31T15:42:07.361-05:00Python: Basics<h2>
2014 Python lecture. Part I </h2>
<br />
The introduction to Python I'm giving at IA-UNAM is accessible here:<br />
<br />
<a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Python.ipynb">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Python.ipynb</a><br />
<br />
I will modify this notebook during the lecture (August 2014), so reload it to have the latest version.<br />
<br />
The topics of this first lecture are:<br />
<ul>
<li>Using python as a calculator</li>
<li>assignments</li>
<li>comments</li>
<li>types</li>
<li>complex numbers</li>
<li>booleans</li>
<li>printing strings</li>
<li>strings</li>
<li>Tuples, lists and dictionaries</li>
<li>Blocks</li>
<li>List and dictionary comprehension</li>
<li>Functions, procedures</li>
<li>Scripting</li>
<li>Importing libraries</li>
</ul>
If you want to have an interactive session with this lecture using the ipython notebook facilities, follow the link above and download the ipynb file (download button at the right top of the web page). Save the file in a directory from where you execute the following (you must have a recent version of ipython installed):<br />
<i>ipython notebook</i><br />
It should open a new tab in your web browser, with the list of ipynb files in the directory. Click on the one you want, will open a new tab similar to the first one, but this one is executed on YOUR computer, it means you are able to interact with the commands. You can change the commands, and execute a cell by SHIFT-ENTER. You can add comments in new cells, and save the result.<br />
<br />
Any comments are welcome.Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-70511151230226588492014-08-07T17:58:00.001-05:002016-05-31T15:42:48.622-05:00Brief introduction to Python<h2>
2014 Python lecture. Part 0 </h2>
<br />
Back to the Python lecture, I want to share here the very quick introduction I gave before starting to play with python: <a href="https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Intro_1.pdf">https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Intro_1.pdf</a><a href="https://drive.google.com/file/d/0B4A0EADFiYFpOFVrdVBHLU5VLUE/edit?usp=sharing"> </a><br />
You may want to install python from Ureka from this site: <a href="http://ssb.stsci.edu/ureka/">http://ssb.stsci.edu/ureka/</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-89486231296631547912014-03-19T22:00:00.001-06:002014-03-19T22:00:36.047-06:00Using ipython Notebook to teach scientific pythonA very good and efficient way to teach python and python related tools, is to use ipython Notebook: <a href="http://ipython.org/notebook.html">http://ipython.org/notebook.html</a><br />
<br />
As example of this use, the following link is a collection of lectures on python, numpy, scipy, matplotlib, use of Fortran from python, etc:<br />
<a href="https://github.com/jrjohansson/scientific-python-lectures">https://github.com/jrjohansson/scientific-python-lectures</a><br />
Enjoy them.<br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-75696156316832485232013-06-28T11:31:00.001-05:002014-08-08T09:27:30.128-05:00Installing python and nice lectures on scientific Python.Since more than one year without any message!... And the new one is almost nothing from me, just links to good pages.<br />
<br />
Two easy ways to install python+ipython+numpy+matplotlib+scipy:<br />
<br />
<b>Anaconda</b> from continuum (but only 64bit version for OSX, which can be a problem for MySQLdb):<br />
Ask for the academic licence if you can.<br />
<a href="http://continuum.io/">http://continuum.io/</a><br />
Once installed, you will need to add the anaconda/bin directory to your PATH and the anaconda directory to your PYTHONPATH.<br />
<br />
<b>Canopy</b> from Entought:<br />
<a href="https://www.enthought.com/products/canopy/">https://www.enthought.com/products/canopy/</a><br />
Once installed, you must setup the virtual environment by adding to your .tcshrc:<br />
setenv VIRTUAL_ENV /Users/YOURNAME/Library/Enthought/Canopy_32bit/User<br />
setenv PATH $VIRTUAL_ENV/bin:$PATH<br />
<br />
!!! Warning !!! When using this virtual environment, don't install using pip with the --user option. Install directly, as for example:<br />
pip install pyfits<br />
<br />
<b><span style="color: red;">UPDATE: Another (better) package comes from STSCI, it's UREKA: <a href="http://ssb.stsci.edu/ureka/">http://ssb.stsci.edu/ureka/</a></span></b><br />
<br />
Here follows good pages to learn python. Those pages are made using Notebook, which is very efficient to show/share python programs.<br />
<br />
Here are the links:<br />
<br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-0-Scientific-Computing-with-Python.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-0-Scientific-Computing-with-Python.ipynb</a><br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-1-Introduction-to-Python-Programming.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-1-Introduction-to-Python-Programming.ipynb</a><br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-2-Numpy.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-2-Numpy.ipynb</a><br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-3-Scipy.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-3-Scipy.ipynb</a><br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-4-Matplotlib.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-4-Matplotlib.ipynb</a><br />
<a href="http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-6B-HPC.ipynb">http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-6B-HPC.ipynb</a><br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-56100476484505572682012-05-05T22:37:00.004-05:002014-09-03T16:22:03.487-05:00Playing with arrays: slicing, sorting, filtering, where function, etc.<h3>
You can also read the more recent post on Numpy here: <a href="http://python-astro.blogspot.mx/2014/08/introduction-to-numpy.html">http://python-astro.blogspot.mx/2014/08/introduction-to-numpy.html</a></h3>
<br />
Most of the data we have to manage are in form of arrays. That's why it's quite always necessary to import the numpy library to easily deal with tables, arrays, images, cubes, etc...<br />
In this post I will discuss some of the methods and functions one may have to use in this context.<br />
<br />
In Python, one can use lists, tuples and dictionaries to put different elements together. They even can contain elements of different types. But nothing better than numpy arrays to really "play" with the data, extract subsets, combine them using arithmetic operations.<br />
<br />
We already discuss a little some array commands in a previous post: <a href="http://python-astro.blogspot.mx/2012/02/read-ascii-file-cont.html">read-ascii-file-cont</a>. Here we are more complete, but not exhaustive, as numpy is a whole world... Have a look at the reference guide: <a href="http://docs.scipy.org/doc/numpy/reference/">http://docs.scipy.org/doc/numpy/reference/</a> and other references at the end of this post.<br />
<h3>
Create arrays</h3>
<div style="background-color: #fce5cd;">
a = np.array([1, 2, 3, 7, 5])</div>
<div style="background-color: #fce5cd;">
a</div>
<div style="background-color: #cfe2f3;">
array([1, 2, 3, 7, 5])<span style="background-color: white;"> </span></div>
<br />
numpy's arrays cannot contain values of different types, the values are transformed if necessary, from integer into real, and from real into string:<br />
<div style="background-color: #fce5cd;">
b = np.array([1, 2, 3, 4.])</div>
<div style="background-color: #fce5cd;">
b</div>
<div style="background-color: #cfe2f3;">
array([ 1., 2., 3., 4.])</div>
<div style="background-color: #fce5cd;">
c = np.array([1, 2, 3., '4'])</div>
<div style="background-color: #fce5cd;">
c</div>
<div style="background-color: #cfe2f3;">
array(['1', '2', '3', '4'], <br />
dtype='|S1')</div>
<br />
shapes and sizes of the arrays are obtained using method and functions:<br />
<div style="background-color: #fce5cd;">
print c.shape</div>
<span style="background-color: #cfe2f3;">(4,)</span><br />
<div style="background-color: #fce5cd;">
print len(c)</div>
<span style="background-color: #cfe2f3;">4</span><br />
<br />
Different ways to create arrays:<br />
<div style="background-color: #fce5cd;">
a = np.arange(1, 10, 0.5)</div>
<div style="background-color: #fce5cd;">
print a</div>
<div style="background-color: #d0e0e3;">
<span style="background-color: #cfe2f3;">[ 1. 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5 6. 6.5 7. 7.5 8.</span><br />
<span style="background-color: #cfe2f3;"> 8.5 9. 9.5]</span></div>
<br />
TAKE CARE, THE LATEST ELEMENT IS NOT WHAT YOU MAY EXPECT...<br />
It's because the first element in Python is indexed by 0. The simple use of np.arange is :<br />
<div style="background-color: #fce5cd;">
a = np.arange(10)</div>
<div style="background-color: #fce5cd;">
print a</div>
<div style="background-color: #d0e0e3;">
[0 1 2 3 4 5 6 7 8 9]</div>
<i>a</i> starts at 0, and has 10 elements, that's ok.<br />
<br />
Easy to create linearly spaced arrays:<br />
<div style="background-color: #fce5cd;">
c = np.linspace(0, 1, 4)</div>
<div style="background-color: #fce5cd;">
print c</div>
<div style="background-color: #cfe2f3;">
[ 0. 0.33333333 0.66666667 1. ]</div>
or log spaced:<br />
<br />
<div style="background-color: #fce5cd;">
c = np.logspace(0, 1, 4)<br />
c</div>
<div style="background-color: #d0e0e3;">
array([ 1. , 2.15443469, 4.64158883, 10. ])</div>
which is actually the same as :<br />
<div style="background-color: #fce5cd;">
c = 10**np.logspace(0, 1, 4)</div>
<br />
To create arrays of 0. or 1.:<br />
<div style="background-color: #fce5cd;">
a = np.zeros(10)</div>
<div style="background-color: #fce5cd;">
b = np.ones((5, 4))</div>
Here <i>b</i> is a 2D array.<br />
<div style="background-color: #fce5cd;">
b.shape</div>
<div style="background-color: #cfe2f3;">
(5, 4)</div>
<br />
You can use a list (or a tuple, or an array) and replicate it:<br />
<div style="background-color: #fce5cd;">
a = np.array([1, 2, 3])</div>
<div style="background-color: #fce5cd;">
b = np.tile(a, 5)</div>
<div style="background-color: #fce5cd;">
print b </div>
<div style="background-color: #cfe2f3;">
[1 2 3 1 2 3 1 2 3 1 2 3 1 2 3]</div>
<div style="background-color: #fce5cd;">
c = np.tile(a, (5, 1))</div>
<div style="background-color: #fce5cd;">
print c</div>
<div style="background-color: #cfe2f3;">
[[1 2 3]<br />
[1 2 3]<br />
[1 2 3]<br />
[1 2 3]<br />
[1 2 3]]</div>
<br />
You can also create 2D arrays from 1D vectors:<br />
<div style="background-color: #fce5cd;">
x = np.linspace(-1, 1, 100)</div>
<div style="background-color: #fce5cd;">
y = x</div>
<div style="background-color: #fce5cd;">
X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel</div>
<div style="background-color: #fce5cd;">
print X.shape </div>
<div style="background-color: #cfe2f3;">
(100, 100)</div>
Outer products are obtain using... np.outer:<br />
<div style="background-color: #fce5cd;">
a = np.array([1, 2, 3])<br />
b = np.array([1, 10, 100])</div>
<div style="background-color: #fce5cd;">
np.outer(a, b)</div>
<span style="background-color: #cfe2f3;">array([[ 1, 10, 100],</span><br />
<span style="background-color: #cfe2f3;"> [ 2, 20, 200],</span><br />
<span style="background-color: #cfe2f3;"> [ 3, 30, 300]])</span><br />
<span style="background-color: #fce5cd;">np.outer(b, a)</span><br />
<div style="background-color: #cfe2f3;">
array([[ 1, 2, 3],<br />
[ 10, 20, 30],<br />
[100, 200, 300]])</div>
<h3>
<b>Broadcasting</b></h3>
Python numpy is able to add dimensions to arrays to perform operations:<br />
<div style="background-color: #fce5cd;">
a = np.array([1, 2, 3, 4])</div>
<div style="background-color: #fce5cd;">
b = np.ones((6, 4))</div>
<div style="background-color: #fce5cd;">
a * b</div>
<div style="background-color: #cfe2f3;">
array([[ 1., 2., 3., 4.],<br />
[ 1., 2., 3., 4.],<br />
[ 1., 2., 3., 4.],<br />
[ 1., 2., 3., 4.],<br />
[ 1., 2., 3., 4.],<br />
[ 1., 2., 3., 4.]])</div>
Nice page on this:<br />
<a href="http://www.scipy.org/EricsBroadcastingDoc">http://www.scipy.org/EricsBroadcastingDoc</a><br />
<h3 style="font-weight: normal;">
<b>Indexing and Slicing</b></h3>
The access to elements of arrays is done using []:<br />
<div style="background-color: #fce5cd;">
a = np.arange(10)</div>
<div style="background-color: #fce5cd;">
a[5]</div>
<div style="background-color: #cfe2f3;">
0 </div>
<div style="background-color: #fce5cd;">
a[-1] # last elements</div>
<span style="background-color: #cfe2f3;">9</span><br />
<br />
In case of N-dims arrays, one can extract slides using :<br />
<div style="background-color: #fce5cd;">
a = np.array([1, 2, 3, 4, 5])</div>
<div style="background-color: #fce5cd;">
b = np.array([1, 10, 100, 1000]) </div>
<div style="background-color: #fce5cd;">
c = np.outer(a, b)</div>
<div style="background-color: #fce5cd;">
c.shape </div>
<div style="background-color: #cfe2f3;">
(5, 4)</div>
<div style="background-color: #fce5cd;">
print c </div>
<div style="background-color: #cfe2f3;">
[[ 1 10 100 1000]<br />
[ 2 20 200 2000]<br />
[ 3 30 300 3000]<br />
[ 4 40 400 4000]<br />
[ 5 50 500 5000]]</div>
<div style="background-color: #fce5cd;">
print c[:,2]</div>
<div style="background-color: #cfe2f3;">
[100 200 300 400 500] </div>
<div style="background-color: #fce5cd;">
print c[-1,:]</div>
<div style="background-color: #cfe2f3;">
[ 5 50 500 5000] </div>
<div style="background-color: #fce5cd;">
print c[-1, -1]</div>
<div style="background-color: #cfe2f3;">
5000</div>
<br />
There is a lot of methods in the array object, the best page to learn more is there: <a href="http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object">http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object</a><br />
<h3>
Filtering</h3>
Sometimes we need to extract elements from arrays following some criteria. There is basically two ways to do this: defining a set of indices where the condition is completed (a la WHERE in IDL), or defining a boolean mask.<br />
<br />
Let's first define a 2D array made of 10 times 1000 random values:<br />
<div style="background-color: #fce5cd;">
a = np.random.random((10, 1000))</div>
<br />
We want to extract the values where the 2nd and the 4th 1000-elements vectors are greater than 0.5.<br />
<br />
Using the <i>numpy.where</i> function:<br />
<div style="background-color: #fce5cd;">
w1 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))</div>
as the result is a tuple of indices, and in this case the dimension of the result is 1, better extract on the fly the array of indices from the tuple:<br />
<div style="background-color: #fce5cd;">
w2 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))[0]</div>
<div style="background-color: #fce5cd;">
len(w2)</div>
<span style="background-color: #cfe2f3;">248 # your result may differ form this value, but must be close to 250.</span><br />
<br />
We can now reduce the initial array to the desired values:<br />
<div style="background-color: #fce5cd;">
b = a[:, w2]</div>
<div style="background-color: #fce5cd;">
b.shape</div>
<div style="background-color: #cfe2f3;">
(10, 248) </div>
<br />
One can also build an array of boolean:<br />
<div style="background-color: #fce5cd;">
mask = (a[1,:] > 0.5) & (a[3,:] > 0.5)</div>
which is actually similar to<br />
<div style="background-color: #fce5cd;">
mask2 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5), True, False)</div>
<br />
It can be used the same way as before:<br />
<div style="background-color: #fce5cd;">
b = a[:, mask]</div>
<div style="background-color: #fce5cd;">
b.shape</div>
<div style="background-color: #cfe2f3;">
(10, 248) </div>
<br />
These latter cases result in a 1000-element array, filled with True and False. To know the number of True values, just sum the array:<br />
<div style="background-color: #fce5cd;">
mask.sum()</div>
<span style="background-color: #cfe2f3;">248</span><span style="background-color: #cfe2f3;"> # your result may differ form this value, but must be close to 250.</span><span style="background-color: #cfe2f3;"> </span><br />
<br />
One big advantage of the mask technique is that you can combine different masks:<br />
<div style="background-color: #fce5cd;">
mask1 = a[1,:] > 0.5</div>
<div style="background-color: #fce5cd;">
mask2 = a[3,:] > 0.5</div>
<div style="background-color: #fce5cd;">
mask_total = mask1 & mask2</div>
<div style="background-color: #fce5cd;">
b = a[:, mask_total]</div>
<h3>
Sorting</h3>
<div style="background-color: #fce5cd;">
a = np.array([1,2,45,2,1,46,7,-1])</div>
<div style="background-color: #fce5cd;">
b1 = np.sort(a)</div>
<div style="background-color: #fce5cd;">
ib = np.argsort(a)</div>
<div style="background-color: #fce5cd;">
b2 = a[ib]</div>
<div style="background-color: #fce5cd;">
print b1</div>
<div style="background-color: #fce5cd;">
<span style="background-color: #cfe2f3;">[-1 1 1 2 2 7 45 46]</span> </div>
<div style="background-color: #fce5cd;">
print b2 </div>
<span style="background-color: #cfe2f3;">[-1 1 1 2 2 7 45 46]</span><br />
<h3>
More on numpy arrays</h3>
<a href="http://scipy.org/Numpy_Example_List_With_Doc">http://scipy.org/Numpy_Example_List_With_Doc</a> <br />
<a href="http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object">http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object</a><br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com7tag:blogger.com,1999:blog-8801264936720939312.post-46196434408544734242012-04-08T16:55:00.000-05:002014-09-11T07:37:51.226-05:00Plotting<br />
Don't miss the update: "2014 Python lecture on Matplotlib" on this blog: <a href="http://python-astro.blogspot.mx/2014/09/2014-python-lecture.html">http://python-astro.blogspot.mx/2014/09/2014-python-lecture.html </a><br />
<br />
Plotting is certainly one of the most common task one would do in Python (at least for astronomers).<br />
There is various libraries to draw plots in Python, but the mostly used and powerful is perhaps <i>matplotlib.</i><br />
The plotting part is pyplot, so before any use, one must import it, with the (widely used) alias <i>plt</i>:<br />
<div style="background-color: #fce5cd;">
import matplotlib.pyplot as plt</div>
<br />
The best to see how it works is to dive into the website where there is a lot of examples of plots with the code used to generate them: <a href="http://matplotlib.sourceforge.net/">http://matplotlib.sourceforge.net/</a>, especially the gallery: <a href="http://matplotlib.sourceforge.net/gallery.html">http://matplotlib.sourceforge.net/gallery.html</a><br />
In the following we will present a few examples to help you to start.<br />
<h3>
X-Y plot</h3>
Let's say you have 2 vectors x and y and wanna plot y vs. x. Here are some examples of how to do this.<br />
First generate x and y (don't forget to <i>import numpy as np</i>):<br />
<div style="background-color: #fce5cd;">
x = np.linspace(0., 4 * np.pi, 100)</div>
<div style="background-color: #fce5cd;">
y = np.sin(x)</div>
<span style="background-color: #fce5cd;">plt.plot(x, y)</span><br />
<span style="background-color: #fce5cd;">plt.show()</span><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVFrtD6YUe5ji3jCObDTQNmVgacxAX89gS2haT3bEfK-VjK_jeHHizn_m2aQjjff42kRl5b_wiG_u4B3U_GjPgM5a7wxal5o9R9GUtbU659-lCgzISDj1GW1TR8qFf7lmWaQQYZ88aD9I/s1600/sin.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVFrtD6YUe5ji3jCObDTQNmVgacxAX89gS2haT3bEfK-VjK_jeHHizn_m2aQjjff42kRl5b_wiG_u4B3U_GjPgM5a7wxal5o9R9GUtbU659-lCgzISDj1GW1TR8qFf7lmWaQQYZ88aD9I/s320/sin.png" height="240" width="320" /></a></div>
<br />
This latest command is not always necessary, if you run ipython with the --pylab option.<br />
As you can see, python automatically define the axis ranges and draw the plot using a default blue color.<br />
If you want to overplot another function, just call plt.plot again:<br />
<div style="background-color: #fce5cd;">
plt.plot(x, y**2)</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjXzgLaXOCJ8_RsNuMJaGCgFmhh2gXssStjlgB_z7YlEbuHnzqD7uaFlZEaKV9KDIFRpBQc6hmCOJEnlfdcutA1cJiZZJC8Q7osXXwDgmWnGM_U2D6LQOPegb57fJEvP8HXMnUjNHWUaYM/s1600/sin2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjXzgLaXOCJ8_RsNuMJaGCgFmhh2gXssStjlgB_z7YlEbuHnzqD7uaFlZEaKV9KDIFRpBQc6hmCOJEnlfdcutA1cJiZZJC8Q7osXXwDgmWnGM_U2D6LQOPegb57fJEvP8HXMnUjNHWUaYM/s320/sin2.png" height="240" width="320" /></a></div>
The figure can be cleaned before another plot using<br />
<div style="background-color: #fce5cd;">
plt.clf()</div>
You may need to draw various figure in different windows at the same time. Every call to<br />
<div style="background-color: #fce5cd;">
plt.figure()</div>
will start a new figure where the next plot(s) will be drawn. The figures (each one in a separate window) are identified with numbers, the first one being 1.<br />
When calling <i>plt.figure(N)</i>, a new figure is created if figure N doesn't exist, and focus will be on figure N otherwise.<br />
<br />
The color of the line is defined using <i>c</i> or <i>color</i> keyword:<br />
<div style="background-color: #fce5cd;">
plt.plot(x, y**3, color = 'red')</div>
<div style="background-color: #fce5cd;">
plt.plot(x, y**4, c='b')</div>
<table border="1" class="docutils"><thead valign="bottom">
<tr><th class="head">Abbreviation</th>
<th class="head">Color</th>
</tr>
</thead>
<tbody valign="top">
<tr><td>b</td>
<td>blue</td>
</tr>
<tr><td>g</td>
<td>green</td>
</tr>
<tr><td>r</td>
<td>red</td>
</tr>
<tr><td>c</td>
<td>cyan</td>
</tr>
<tr><td>m</td>
<td>magenta</td>
</tr>
<tr><td>y</td>
<td>yellow</td>
</tr>
<tr><td>k</td>
<td>black</td>
</tr>
<tr><td>w</td>
<td>white</td></tr>
</tbody></table>
<br />
<br />
Symbols and line style can also be easily defined:<br />
<div style="background-color: #fce5cd;">
plt.plot(x, y**3, color='r', marker='o', linestyle=':')<span style="background-color: white;"></span> </div>
<table border="1" class="docutils"><thead valign="bottom">
<tr><th class="head">Symbol</th>
<th class="head">Description</th>
</tr>
</thead>
<tbody valign="top">
<tr><td><tt class="docutils literal"><span class="pre">-</span></tt></td>
<td>solid line</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">--</span></tt></td>
<td>dashed line</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">-.</span></tt></td>
<td>dash-dot line</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">:</span></tt></td>
<td>dotted line</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">.</span></tt></td>
<td>points</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">,</span></tt></td>
<td>pixels</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">o</span></tt></td>
<td>circle symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">^</span></tt></td>
<td>triangle up symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">v</span></tt></td>
<td>triangle down symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre"><</span></tt></td>
<td>triangle left symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">></span></tt></td>
<td>triangle right symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">s</span></tt></td>
<td>square symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">+</span></tt></td>
<td>plus symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">x</span></tt></td>
<td>cross symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">D</span></tt></td>
<td>diamond symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">d</span></tt></td>
<td>thin diamond symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">1</span></tt></td>
<td>tripod down symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">2</span></tt></td>
<td>tripod up symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">3</span></tt></td>
<td>tripod left symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">4</span></tt></td>
<td>tripod right symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">h</span></tt></td>
<td>hexagon symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">H</span></tt></td>
<td>rotated hexagon symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">p</span></tt></td>
<td>pentagon symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">|</span></tt></td>
<td>vertical line symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">_</span></tt></td>
<td>horizontal line symbols</td>
</tr>
<tr><td><tt class="docutils literal"><span class="pre">steps</span></tt></td>
<td>use gnuplot style ‘steps’ # kwarg only</td></tr>
</tbody></table>
<br />
Other useful line properties:<br />
<br />
<table border="1" class="docutils"><thead valign="bottom">
<tr><th class="head">Property</th>
<th class="head">Value</th>
</tr>
</thead>
<tbody valign="top">
<tr><td>alpha</td>
<td>alpha transparency on 0-1 scale</td>
</tr>
<tr><td>antialiased</td>
<td>True or False - use antialised rendering</td>
</tr>
<tr><td>color</td>
<td>matplotlib color arg</td>
</tr>
<tr><td>data_clipping</td>
<td>whether to use numeric to clip data</td>
</tr>
<tr><td>label</td>
<td>string optionally used for legend</td>
</tr>
<tr><td>linestyle</td>
<td>one of <tt class="docutils literal"><span class="pre">-</span></tt> <tt class="docutils literal"><span class="pre">:</span></tt> <tt class="docutils literal"><span class="pre">-.</span></tt> <tt class="docutils literal"><span class="pre">-</span></tt></td>
</tr>
<tr><td>linewidth</td>
<td>float, the line width in points</td>
</tr>
<tr><td>marker</td>
<td>one of <tt class="docutils literal"><span class="pre">+</span></tt> <tt class="docutils literal"><span class="pre">,</span></tt> <tt class="docutils literal"><span class="pre">o</span></tt> <tt class="docutils literal"><span class="pre">.</span></tt> <tt class="docutils literal"><span class="pre">s</span></tt> <tt class="docutils literal"><span class="pre">v</span></tt> <tt class="docutils literal"><span class="pre">x</span></tt> <tt class="docutils literal"><span class="pre">></span></tt> <tt class="docutils literal"><span class="pre"><</span></tt>, etc</td>
</tr>
<tr><td>markeredgewidth</td>
<td>line width around the marker symbol</td>
</tr>
<tr><td>markeredgecolor</td>
<td>edge color if a marker is used</td>
</tr>
<tr><td>markerfacecolor</td>
<td>face color if a marker is used</td>
</tr>
<tr><td>markersize</td>
<td>size of the marker in points</td>
</tr>
</tbody>
</table>
<br />
If you need more control on the markers, better use scatter. For example, if one need ot change the size of the symbol according to a function:<br />
<div style="background-color: #fce5cd;">
plt.scatter(x, y**3, c=abs(y), marker='s', s=10+abs(y)*100, edgecolors='none')</div>
<br />
The labels are set after the plot is done:<br />
<div style="background-color: #fce5cd;">
plt.xlabel('X')</div>
<div style="background-color: #fce5cd;">
plt.ylabel('Some trig function')</div>
LaTex fans are welcome:<br />
<div style="background-color: #fce5cd;">
plt.title(r'Example of LaTex $\alpha_{\beta}$') #Notice the r before the string</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj7ODLeoQCBXYLb6VX4ovWLOCJ7LVwRi0vneeZ82zNLZhzJiul_V3I0Wb7uHGtVQQNElITtDjp7ImUxa7eMR9peRhxhSyRvOa97843mHR7z1RcAB1GfmG-R6d9bmDWGBKmAo6jTN4aV42A/s1600/sin3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj7ODLeoQCBXYLb6VX4ovWLOCJ7LVwRi0vneeZ82zNLZhzJiul_V3I0Wb7uHGtVQQNElITtDjp7ImUxa7eMR9peRhxhSyRvOa97843mHR7z1RcAB1GfmG-R6d9bmDWGBKmAo6jTN4aV42A/s320/sin3.png" height="240" width="320" /></a></div>
<br />
One can change the axis ranges afterward:<br />
<div style="background-color: #fce5cd;">
plt.xlim((0, 10))</div>
<div style="background-color: #fce5cd;">
plt.ylim((-2, 2))</div>
<h3>
Multiple plots</h3>
To draw multiple plots of the same figure:<br />
<div style="background-color: #fce5cd;">
plt.subplot(N_y, N_x, N)</div>
<br />
<div style="background-color: #fce5cd;">
for i in np.arange(9):</div>
<div style="background-color: #fce5cd;">
plt.subplot(3,3,i+1)</div>
<span style="background-color: #fce5cd;"> plt.plot(x, y**i)</span><br />
<span style="background-color: #fce5cd;"> plt.ylim((-2, 2))</span><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgE9Sat9b1cg-oMXs2337n631ta7BF1Q1H6Fyvh3VZRAGAERDQly44JUBnNBP5lkjfoEH9s54UV5StDDQCKkDO2foc3JIGTXTZG6TkOv6XDy2gUA6SpkBpOVn5wIXGTW72naSaMQ57TNg8/s1600/sin4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgE9Sat9b1cg-oMXs2337n631ta7BF1Q1H6Fyvh3VZRAGAERDQly44JUBnNBP5lkjfoEH9s54UV5StDDQCKkDO2foc3JIGTXTZG6TkOv6XDy2gUA6SpkBpOVn5wIXGTW72naSaMQ57TNg8/s320/sin4.png" height="238" width="320" /></a></div>
<h3>
log plots</h3>
using <i>plt.semilogx, plt.semilogy</i> and <i>plt.loglog</i> <br />
<div style="background-color: #fce5cd;">
plt.loglog(x)</div>
<div style="background-color: #fce5cd;">
plt.grid(True, which='minor')</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglTfAMBpxKaOmmvY_YgjYByQEoRHqfMYKixHZfOcIeTYjpqQHkn7UmvK5b2PIhyphenhyphentaK4i3lXE-_YAvk8yQMI8t-DOl3v6hzh0Fy2t7_Wc_yS-xAmdD87vhKhwKBlq5C3szywXKWGSJtM0A/s1600/log.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglTfAMBpxKaOmmvY_YgjYByQEoRHqfMYKixHZfOcIeTYjpqQHkn7UmvK5b2PIhyphenhyphentaK4i3lXE-_YAvk8yQMI8t-DOl3v6hzh0Fy2t7_Wc_yS-xAmdD87vhKhwKBlq5C3szywXKWGSJtM0A/s320/log.png" height="240" width="320" /></a></div>
<h3>
Contours </h3>
Contour plots are done with <i>plt.contour</i> and <i>plt.contourf</i>.<br />
<div style="background-color: #fce5cd;">
x = np.linspace(-1, 1, 100)</div>
<div style="background-color: #fce5cd;">
y = x</div>
<div style="background-color: #fce5cd;">
X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel</div>
<div style="background-color: #fce5cd;">
dist = (X**2 + Y**2)**0.5</div>
<div style="background-color: #fce5cd;">
plt.contourf(X, Y, dist)</div>
<div style="background-color: #fce5cd;">
plt.colorbar()</div>
<div style="background-color: #fce5cd;">
CS = plt.contour(X, Y, dist, colors ='black', linewidths = 4)</div>
<span style="background-color: #fce5cd;">plt.clabel(CS) #to print label on each contour</span><br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifBXz0NmUxtq4mZgzF_sRq0yuqGnyPoWsOFbHSB75Czi4ikEWy9wRvce1TEgp2bskWFUFpt3jjJek_WV06WulZ-QRVe5kYpTcMMDbg0KGfiCt_Yegd1EWLwqWFHP3cdDt1fKk9N-X1zW8/s1600/contour1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifBXz0NmUxtq4mZgzF_sRq0yuqGnyPoWsOFbHSB75Czi4ikEWy9wRvce1TEgp2bskWFUFpt3jjJek_WV06WulZ-QRVe5kYpTcMMDbg0KGfiCt_Yegd1EWLwqWFHP3cdDt1fKk9N-X1zW8/s320/contour1.png" height="240" width="320" /></a></div>
<br />
<h3>
Saving the plot</h3>
The result of the plot can be save in PDF, EPS, JPG, BMP format, using:<br />
<span style="background-color: #fce5cd;">plt.savefig('fig1.pdf')</span><br />
<h3>
More on plotting:</h3>
<a href="http://scipy-lectures.github.com/intro/matplotlib/matplotlib.html">http://scipy-lectures.github.com/intro/matplotlib/matplotlib.html</a><br />
<br />
<br />Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com1tag:blogger.com,1999:blog-8801264936720939312.post-19328347708608841682012-03-25T11:21:00.001-06:002014-09-03T16:27:14.495-05:00Play with FITS files<h3>
The same lecture is now in the Notebook format, see there: <a href="http://python-astro.blogspot.mx/2014/09/interacting-with-files-reading-writing.html">http://python-astro.blogspot.mx/2014/09/interacting-with-files-reading-writing.html</a></h3>
<h3>
<br /></h3>
<h3>
What is the FITS format?</h3>
The FITS format is the most popular way to save and interchange astronomical data. The files are organized in units each of which contains a human readable header and a data. This structure is refereed as HDUs (Header/DATA Unit). <br />
<br />
A FITS file can contain one or more HDUs, the first of which is called "primary" and the rest are called "extensions". The primary HDU usually contains 1D spectrum, 2D image or 3D data cube, although any dimension from 0 to 999 are possible. The data are 1, 2 or 4 bytes integers or 4 or 8 bytes real numbers. <br />
<br />
The extensions can contain or arrays as in the primary HDU or ascii tables or binary tables. <br />
If a FITS file contains only tables, it primary HDU does not contain data, but only header. <br />
<br />
Both headers and data in a FITS file are organized in blocs of 2880 bytes. The header contain 80 bytes lines each of which consists of a keyword of 8 bytes followed in most of the cases by '= ' in the position 9 and 10 and then the value of the keyword. The rest of the line is a comment string beginning with '/'. Each header begins with the following lines<br />
<br />
SIMPLE = T / file conforms to FITS standard<br />
BITPIX = 16 / number of bits per data pixel<br />
NAXIS = 2 / number of data axes<br />
NAXIS1 = 440 / length of data axis 1<br />
NAXIS2 = 300 / length of data axis 2<br />
<br />
which defines the format of the file as standard FITS, the data format and the dimensions of the stored data. <br />
<br />
One block of 2880 bytes contains 36 lines of 80 characters per line. The header can have several blocks of 36 lines. The last block is identified by the presence of the keyword 'END' The next 2880 bytes block contains the first part of the data. The empty lines after 'END' keyword are filled with blanks and the unused bytes from the end of the data to the end of the 2880 bytes block are filled with NULLs.<br />
<br />
Full description of the FITS format can be found at <a href="http://fits.gsfc.nasa.gov/fits_primer.html">http://fits.gsfc.nasa.gov/fits_primer.html</a><br />
<h3>
PyFITS</h3>
As one can see, reading FITS files is not as simple as reading columns of data. Fortunately, there is a python package called PyFITS which makes the task much easier. <br />
<br />
The package can be downloaded from <br />
<br />
<a href="http://www.stsci.edu/institute/software_hardware/pyfits">http://www.stsci.edu/institute/software_hardware/pyfits</a><br />
<br />
The documentation can be found on the same web page. <br />
<br />
Although it is possible to download and install the package it is much easier to install it with <br />
<br />
easy_install pyfits<br />
or<br />
pip pyfits <br />
<br />
Information on easy_install: <a href="http://packages.python.org/distribute/easy_install.html">HERE</a> and on pip: <a href="http://www.pip-installer.org/en/latest/index.html">HERE</a>.<br />
I did that on both Mac and Linux machines and it working nicely. <br />
<h3>
How to read FITS file</h3>
First, as usual, one has to import the pyfits package<br />
<br />
<span style="background-color: #fce5cd;">import pyfits</span><br />
<br />
then the file is read with <br />
<br />
<span style="background-color: #fce5cd;">hdulist = pyfits.open('n10017o.fits')</span><br />
<br />
where I used one of my FITS files from San Pedro Martir echelle spectrograph. The file can be downloaded from <a href="https://docs.google.com/open?id=0B4A0EADFiYFpVXo3QWdhbk5UM09Ld3BvSWl4d2tBUQ">HERE</a>.<br />
<br />
The result hdulist is a list of HDU objects. In the case of a simple file, there is only one primary HDU so the list contains only one element<br />
<br />
<span style="background-color: #fce5cd;">len(hdulist)</span><br />
<span style="background-color: #cfe2f3;">1</span><br />
<br />
The information on what the file contains can be obtained by calling the info() method:<br />
<br />
<span style="background-color: #fce5cd;">hdulist.info()</span><br />
<span style="background-color: #cfe2f3;">Filename: n10017o.fits<br />No. Name Type Cards Dimensions Format<br />0 PRIMARY PrimaryHDU 62 (2154, 2048) int16</span><br />
<br />
The table said that there is only a primary HDU which contains 2154 X 2048 image with data stored in 2 bytes (16 bits) integers.<br />
<br />
As described above, the HDU (header/data unit) contains header and data. The header is a dictionary. To see what keywords were used in the header one can do:<br />
<br />
<span style="background-color: #fce5cd;">hdulist[0].header.keys()</span><br />
<span style="background-color: #cfe2f3;">['SIMPLE',</span><br />
<span style="background-color: #cfe2f3;"> 'BITPIX',</span><br />
<span style="background-color: #cfe2f3;"> 'NAXIS',</span><br />
<span style="background-color: #cfe2f3;"> 'NAXIS1',</span><br />
<span style="background-color: #cfe2f3;"> 'NAXIS2',</span><br />
<span style="background-color: #cfe2f3;"> 'EXTEND',</span><br />
<span style="background-color: #cfe2f3;"> 'COMMENT',</span><br />
<span style="background-color: #cfe2f3;"> 'BZERO',</span><br />
<span style="background-color: #cfe2f3;"> 'BSCALE',</span><br />
<span style="background-color: #cfe2f3;"> 'EXPTIME',</span><br />
<span style="background-color: #cfe2f3;"> 'DETECTOR',</span><br />
<span style="background-color: #cfe2f3;"> 'ORIGIN',</span><br />
<span style="background-color: #cfe2f3;"> 'OBSERVAT',</span><br />
<span style="background-color: #cfe2f3;"> 'TELESCOP',</span><br />
<span style="background-color: #cfe2f3;"> 'LATITUDE',</span><br />
<span style="background-color: #cfe2f3;"> 'LONGITUD',</span><br />
<span style="background-color: #cfe2f3;"> 'ALTITUD',</span><br />
<span style="background-color: #cfe2f3;"> 'SECONDAR',</span><br />
<span style="background-color: #cfe2f3;"> 'TIMEZONE',</span><br />
<span style="background-color: #cfe2f3;"> 'OBSERVER',</span><br />
<span style="background-color: #cfe2f3;"> 'OBJECT',</span><br />
<span style="background-color: #cfe2f3;"> 'INSTRUME',</span><br />
<span style="background-color: #cfe2f3;"> 'GAINMODE',</span><br />
<span style="background-color: #cfe2f3;"> 'FILTER',</span><br />
<span style="background-color: #cfe2f3;"> 'IMGTYPE',</span><br />
<span style="background-color: #cfe2f3;"> 'EQUINOX',</span><br />
<span style="background-color: #cfe2f3;"> 'ST',</span><br />
<span style="background-color: #cfe2f3;"> 'UT',</span><br />
<span style="background-color: #cfe2f3;"> 'JD',</span><br />
<span style="background-color: #cfe2f3;"> 'DATE-OBS',</span><br />
<span style="background-color: #cfe2f3;"> 'CCDSUM',</span><br />
<span style="background-color: #cfe2f3;"> 'RA',</span><br />
<span style="background-color: #cfe2f3;"> 'DEC',</span><br />
<span style="background-color: #cfe2f3;"> 'AH',</span><br />
<span style="background-color: #cfe2f3;"> 'AIRMASS',</span><br />
<span style="background-color: #cfe2f3;"> 'TMMIRROR',</span><br />
<span style="background-color: #cfe2f3;"> 'TSMIRROR',</span><br />
<span style="background-color: #cfe2f3;"> 'TAIR',</span><br />
<span style="background-color: #cfe2f3;"> 'XTEMP',</span><br />
<span style="background-color: #cfe2f3;"> 'HUMIDITY',</span><br />
<span style="background-color: #cfe2f3;"> 'ATMOSBAR',</span><br />
<span style="background-color: #cfe2f3;"> 'WIND',</span><br />
<span style="background-color: #cfe2f3;"> 'WDATE',</span><br />
<span style="background-color: #cfe2f3;"> 'DATE',</span><br />
<span style="background-color: #cfe2f3;"> 'NAMPS',</span><br />
<span style="background-color: #cfe2f3;"> 'CCDNAMPS',</span><br />
<span style="background-color: #cfe2f3;"> 'AMPNAME',</span><br />
<span style="background-color: #cfe2f3;"> 'CREATOR',</span><br />
<span style="background-color: #cfe2f3;"> 'VERSION',</span><br />
<span style="background-color: #cfe2f3;"> 'HISTORY']</span><br />
<br />
and to get the value of a given keyword <br />
<br />
<span style="background-color: #fce5cd;">hdulist[0].header['OBJECT']</span><span style="background-color: #cfe2f3;"> </span><br />
<span style="background-color: #cfe2f3;">'BD +59 363'</span><br />
<br />
The header can be printed as it apears in the file by<br />
<br />
<span style="background-color: #fce5cd;">print hdulist[0].header.ascardlist()</span><br />
<span style="background-color: #cfe2f3;">SIMPLE = T / conforms to FITS standard </span><br />
<span style="background-color: #cfe2f3;">BITPIX = 16 / array data type </span><br />
<span style="background-color: #cfe2f3;">NAXIS = 2 / number of array dimensions </span><br />
<span style="background-color: #cfe2f3;">NAXIS1 = 2154 / length of data axis 1 </span><br />
<span style="background-color: #cfe2f3;">NAXIS2 = 2048 / length of data axis 2 </span><br />
<span style="background-color: #cfe2f3;">EXTEND = T </span><br />
<span style="background-color: #cfe2f3;">COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy </span><br />
<span style="background-color: #cfe2f3;">COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H </span><br />
<span style="background-color: #cfe2f3;">BZERO = 32768 / BZERO </span><br />
<span style="background-color: #cfe2f3;">BSCALE = 1 / BSCALE </span><br />
<span style="background-color: #cfe2f3;">EXPTIME = 0.0 / Integration Time, sec. </span><br />
<span style="background-color: #cfe2f3;">DETECTOR= 'e2vm2 E2V-4240' / CCD Type </span><br />
<span style="background-color: #cfe2f3;">ORIGIN = 'UNAM ' / OAN SPM, IA-UNAM </span><br />
<span style="background-color: #cfe2f3;">OBSERVAT= 'SPM ' / Observatory </span><br />
<span style="background-color: #cfe2f3;">TELESCOP= '2.12m ' / Telescope </span><br />
<span style="background-color: #cfe2f3;">LATITUDE= '+31:02:39' / Latitude </span><br />
<span style="background-color: #cfe2f3;">LONGITUD= '115:27:49' / Longitud </span><br />
<span style="background-color: #cfe2f3;">ALTITUD = 2800 / altitud </span><br />
<span style="background-color: #cfe2f3;">SECONDAR= -1 / F/ Secondary type </span><br />
<span style="background-color: #cfe2f3;">TIMEZONE= 8 / Time Zone </span><br />
<span style="background-color: #cfe2f3;">OBSERVER= 'Leonid ' / Observer's Name </span><br />
<span style="background-color: #cfe2f3;">OBJECT = 'BD +59 363' / Object </span><br />
<span style="background-color: #cfe2f3;">INSTRUME= 'Echelle ' / Instrument </span><br />
<span style="background-color: #cfe2f3;">GAINMODE= 1 / Gain factor in the CCD </span><br />
<span style="background-color: #cfe2f3;">FILTER = 'None ' / Filter </span><br />
<span style="background-color: #cfe2f3;">IMGTYPE = 'zero ' / Image Type </span><br />
<span style="background-color: #cfe2f3;">EQUINOX = 2011.7 / Equinox </span><br />
<span style="background-color: #cfe2f3;">ST = '02:19:51.2' / Sideral Time </span><br />
<span style="background-color: #cfe2f3;">UT = '11:28:28' / Universal Time </span><br />
<span style="background-color: #cfe2f3;">JD = 2455803.5 / Julian Date </span><br />
<span style="background-color: #cfe2f3;">DATE-OBS= '2011-08-30' / Observation Date UTM </span><br />
<span style="background-color: #cfe2f3;">CCDSUM = '1 1 ' / Binning [ Cols:Rows ] </span><br />
<span style="background-color: #cfe2f3;">RA = ' 02:13:22.2' / Right Ascension </span><br />
<span style="background-color: #cfe2f3;">DEC = ' 52''14''44.0' / Declination </span><br />
<span style="background-color: #cfe2f3;">AH = ' 00:06:28.0' / Hour Angle </span><br />
<span style="background-color: #cfe2f3;">AIRMASS = 1.073 / Airmass </span><br />
<span style="background-color: #cfe2f3;">TMMIRROR= 0 / Primary Mirror Temperature (celsius degree) </span><br />
<span style="background-color: #cfe2f3;">TSMIRROR= 0 / Secundary Mirror Temperature (celsius degree) </span><br />
<span style="background-color: #cfe2f3;">TAIR = 0 / Internal Telescope Air Temperature (celsius deg</span><br />
<span style="background-color: #cfe2f3;">XTEMP = 14.6 / Exterior Temperature (celsius degree) </span><br />
<span style="background-color: #cfe2f3;">HUMIDITY= 44.0 / % external Humidity </span><br />
<span style="background-color: #cfe2f3;">ATMOSBAR= 731.7 / Atmosferic Presure in mb </span><br />
<span style="background-color: #cfe2f3;">WIND = 'S at 24.1 km/h' / Wind Direction </span><br />
<span style="background-color: #cfe2f3;">WDATE = '11:28:10, 08/30/11' / Weather Acquisition Date (Local time) </span><br />
<span style="background-color: #cfe2f3;">DATE = '2011-08-30T11:28:29' / file creation date (YYYY-MM-DDThh:mm:ss UT) </span><br />
<span style="background-color: #cfe2f3;">NAMPS = 1 / Number of Amplifiers </span><br />
<span style="background-color: #cfe2f3;">CCDNAMPS= 1 / Number of amplifiers used </span><br />
<span style="background-color: #cfe2f3;">AMPNAME = '1 Channel' / Amplifier name </span><br />
<span style="background-color: #cfe2f3;">CREATOR = 'Python Oan ccds' / Name of the software task that created the file</span><br />
<span style="background-color: #cfe2f3;">VERSION = '4.12D ' / Application Software Version </span><br />
<span style="background-color: #cfe2f3;">COMMENT Visit our weather site http://www.astrossp.unam.mx/weather15 </span><br />
<span style="background-color: #cfe2f3;">COMMENT for complete meteorological data of your observation night </span><br />
<span style="background-color: #cfe2f3;">HISTORY bin2fits V1.0 </span><br />
<span style="background-color: #cfe2f3;">HISTORY Programmer: Enrique Colorado [ colorado@astrosen.unam.mx ] </span><br />
<span style="background-color: #cfe2f3;">HISTORY Observatorio Astronomico Nacional -UNAM </span><br />
<span style="background-color: #cfe2f3;">HISTORY V1.00 By Arturo Nunez and Colorado >Ported to Python using pyfits </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.50 By E. Colorado >Added interior mirrors temperatures </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.49 By E. Colorado >Added BIASSEC parameter </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.48 By E. Colorado >Aditional info for autofocus calculations </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.4 By E. Colorado >Now we include timezone, and remove lat. sign </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.3 By E. Colorado >Now we include weather data </span><br />
<span style="background-color: #cfe2f3;">HISTORY V0.2 By E. Colorado >General OAN Working Release </span><br />
<br />
The data in the file are accessible with<br />
<br />
<span style="background-color: #fce5cd;">data = hdulist[0].data</span><br />
<br />
and can be seen with [don't forget to <i>import matplotlib.pyplot as plt</i> before running this]:<br />
<br />
<span style="background-color: #fce5cd;">plt.imshow(data)</span><br />
<br />
A column from the data can be plotted with <br />
<br />
<span style="background-color: #fce5cd;">plt.plot(data[:,1000])</span><br />
<br />
where I am plotting the column number 1000. In the same way a line from the data is plotted with: <br />
<br />
<span style="background-color: #fce5cd;">plt.plot(data[1000,:])</span><br />
<br />
The data are numpy object so all manipulations are available. <br />
<h3>
Some more on displaying images</h3>
The main matplotlib function to display images is imshow. It has several parameters of which te most important are the one which controls the color scheme and the ones which control the dynamic range. <br />
<br />
<span style="background-color: #fce5cd;">plt.imshow(data, cmap=cm.gray, vmin=1000, vmax=10000)</span><br />
<br />
<i>cmap</i> controls the pseudocolor map (same as color= in IDL). The predefined color maps can be seen on <a href="http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html">http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html</a><br />
<br />
<i>vmin</i> and <i>vmax</i> controls the range of values which are mapped in [0,1] range where 0 is black (the darkest color) and 1 is white (the lightest color).<br />
<br />
<i>origin = 'upper' | 'lower'</i>, place the pixel with coordinates 0,0 at the upper or lower corner of the plot<br />
<br />
<i>extent = (xmin,xmax,ymin,ymax)</i> - This parameter defines the numbers written on the axes. No changes on the image. <br />
<h3>
Using FITS tables</h3>
For this example I'll use a spectrum obtain with the high dispersion camera on board of IUE. <br />
The file is opened as usual:<br />
<span style="background-color: #fce5cd;">hdulist = pyfits.open('swp04345.mxhi')</span><br />
<br />
Download the file from <a href="https://docs.google.com/open?id=0B4A0EADFiYFpTlN5Wkd1cU9SVXFNdzg0WDlWV196UQ">THERE</a>.<br />
<br />
but now hdulist has 2 elements (2 header/data units):<br />
<br />
<span style="background-color: #fce5cd;">len(hdulist)</span><br />
<span style="background-color: #cfe2f3;">2</span><br />
<br />
We can see that the primary header has dimension (), son does not contain any data. The data are in the extension.<br />
<br />
<span style="background-color: #fce5cd;">hdulist.info()</span><br />
<span style="background-color: #cfe2f3;">Filename: swp04345.mxhi</span><br />
<span style="background-color: #cfe2f3;">No. Name Type Cards Dimensions Format</span><br />
<span style="background-color: #cfe2f3;">0 PRIMARY PrimaryHDU 421 () uint8</span><br />
<span style="background-color: #cfe2f3;">1 MEHI BinTableHDU 61 60R x 17C [1B, 1I, 1D, 1I, 1D, 1E, 1E, 768E, 768E, 768E, 768I, 768E, 768E, 1I, 1I, 1E, 7E]</span><br />
<br />
The first header contains the minimal infirmation:<br />
<br />
<span style="background-color: #fce5cd;">print hdulist[0].header.ascardlist()[:5]</span><br />
<span style="background-color: #cfe2f3;">SIMPLE = T / Standard FITS Format </span><br />
<span style="background-color: #cfe2f3;">BITPIX = 8 / Binary data </span><br />
<span style="background-color: #cfe2f3;">NAXIS = 0 / Two-dimensional image </span><br />
<span style="background-color: #cfe2f3;">EXTEND = T / Extensions are present </span><br />
<span style="background-color: #cfe2f3;">TELESCOP= 'IUE ' / International Ultraviolet Explorer </span><br />
<br />
The number of axis is 0 which means there is no data block in the primary HDU. The header of the second HDU begins with the keyword XTENSION and with the specification of the data<br />
<br />
<span style="background-color: #fce5cd;">print hdulist[1].header.ascardlist()[:5]</span><br />
<span style="background-color: #cfe2f3;">XTENSION= 'BINTABLE' / Binary table extension </span><br />
<span style="background-color: #cfe2f3;">BITPIX = 8 / Binary data </span><br />
<span style="background-color: #cfe2f3;">NAXIS = 2 / Two-dimensional table array </span><br />
<span style="background-color: #cfe2f3;">NAXIS1 = 16961 / Width of row in bytes </span><br />
<span style="background-color: #cfe2f3;">NAXIS2 = 60 / Number of orders </span> <br />
<br />
To progress further we need to know what is in the table. As usual, the columns have names and type of the stored data. These information can be obtained using the column attribute of hdulist:<br />
<br />
<span style="background-color: #fce5cd;">cols = hdulist[1].columns</span><br />
<br />
<span style="background-color: #fce5cd;">cols.info</span><br />
<span style="background-color: #cfe2f3;"><bound method ColDefs.info of ColDefs(<br /> name = 'ORDER'; format = '1B'; unit = ' '<br /> name = 'NPOINTS'; format = '1I'; unit = ' '<br /> name = 'WAVELENGTH'; format = '1D'; unit = 'ANGSTROM'<br /> name = 'STARTPIX'; format = '1I'; unit = 'PIXEL'<br /> name = 'DELTAW'; format = '1D'; unit = 'ANGSTROM'<br /> name = 'SLIT HEIGHT'; format = '1E'; unit = 'PIXEL'<br /> name = 'LINE_FOUND'; format = '1E'; unit = 'PIXEL'<br /> name = 'NET'; format = '768E'; unit = 'FN'<br /> name = 'BACKGROUND'; format = '768E'; unit = 'FN'<br /> name = 'NOISE'; format = '768E'; unit = 'FN'<br /> name = 'QUALITY'; format = '768I'; unit = ' '<br /> name = 'RIPPLE'; format = '768E'; unit = 'FN'<br /> name = 'ABS_CAL'; format = '768E'; unit = 'ERG/CM2/S/A'<br /> name = 'START-BKG'; format = '1I'; unit = 'PIXEL'<br /> name = 'END-BKG'; format = '1I'; unit = 'PIXEL'<br /> name = 'SCALE_BKG'; format = '1E'; unit = ' '<br /> name = 'COEFF'; format = '7E'; unit = ' '<br />)></span><br />
<br />
the <i>cols.info</i> returns the names of the columns and the information of their format and units.<br />
<br />
The data are available using (this example is NOT the right way of plotting the data, it's just an example) [and don't forget to <i>import numpy as np</i> to have np.arange working]:<br />
<br />
<div style="background-color: #fce5cd;">
data1 = hdulist[1].data</div>
<div style="background-color: #fce5cd;">
DTs = data1.ABS_CAL</div>
<div style="background-color: #fce5cd;">
WLs = data1.WAVELENGTH</div>
<div style="background-color: #fce5cd;">
DWs = data1.DELTAW </div>
<div style="background-color: #fce5cd;">
for WL, DW, DT in zip(WLs, DWs, DTs):</div>
<div style="background-color: #fce5cd;">
plot(WL + np.arange(len(DT)) * DW, DT)</div>
<h3>
Writing FITS files.</h3>
The creation of a FITS file pass through 4 steps. <br />
<br />
1) Creation of numpy array with the data. <br />
<br />
<span style="background-color: #fce5cd;">x = np.arange(100)</span><br />
<br />
2) Creation of the HDU from the data. <br />
<br />
<span style="background-color: #fce5cd;">hdu = pyfits.PrimaryHDU(x)</span><br />
<br />
thus created, the hdu has its basic header and the data. <br />
<br />
<span style="background-color: #fce5cd;">print hdu.header.ascardlist()</span><br />
<span style="background-color: #cfe2f3;">SIMPLE = T / conforms to FITS standard </span><br />
<span style="background-color: #cfe2f3;">BITPIX = 64 / array data type </span><br />
<span style="background-color: #cfe2f3;">NAXIS = 1 / number of array dimensions </span><br />
<span style="background-color: #cfe2f3;">NAXIS1 = 100 </span><br />
<span style="background-color: #cfe2f3;">EXTEND = T </span> <br />
<br />
3) Adding additional keywords to the header. The automatically created header contains only the required minimum of keywords. If additional keywords are needed they are added with<br />
<br />
<span style="background-color: #fce5cd;">hdu.header.update('testkey',0.001,'some test value')</span><br />
<br />
<span style="background-color: #fce5cd;">print hdu.header.ascardlist()</span><br />
<span style="background-color: #cfe2f3;">SIMPLE = T / conforms to FITS standard </span><br />
<span style="background-color: #cfe2f3;">BITPIX = 64 / array data type </span><br />
<span style="background-color: #cfe2f3;">NAXIS = 1 / number of array dimensions </span><br />
<span style="background-color: #cfe2f3;">NAXIS1 = 100 </span><br />
<span style="background-color: #cfe2f3;">EXTEND = T </span><br />
<span style="background-color: #cfe2f3;">TESTKEY = 0.001 / some test value </span> <br />
<br />
4) Once all the keywords are ready, the final HDU list have to be created and written to the file:<br />
<br />
<span style="background-color: #fce5cd;">hdulist = pyfits.HDUList([hdu])</span><br />
<span style="background-color: #fce5cd;">hdulist.writeto('new.fits')</span><br />
<span style="background-color: #fce5cd;">hdulist.close()</span><br />
<h3>
Alternative way </h3>
Another way to deal with FITS tables is to use the ATpy library, look there:<br />
<a href="http://atpy.github.com/index.html">http://atpy.github.com/index.html</a><br />
<h3>
References</h3>
<br />
To learn more on PyFITS: <a href="http://packages.python.org/pyfits/">http://packages.python.org/pyfits/</a>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com2tag:blogger.com,1999:blog-8801264936720939312.post-8354407202885133452012-02-27T12:14:00.000-06:002012-03-12T11:00:54.669-06:00Read an ascii file (cont')<br />
<div style="color: blue;">
<b>In this message we'll see other ways to do the same job than in the previous post: read an ascii file.</b></div>
<br />
In the previous post we read a file of the form:<br />
<br />
0.000000 0.000000 <br />
0.095200 0.095056 <br />
0.190400 0.189251 <br />
0.285599 0.281733 <br />
0.380799 0.371662 <br />
0.475999 0.458227 <br />
0.571199 0.540641 <br />
0.666398 0.618159 <br />
0.761598 0.690079 <br />
0.856798 0.755750<br />
...<br />
using 2 lists, one for each column. It can be easier to manage only one 2D variable containing the whole table. Easier to apply filter on it, for example.<br />
We first read the whole file into a single variable:<br />
<div style="background-color: #fce5cd;">
f2 = open('datas.dat', 'r')<br />
lines = f2.readlines()<br />
f2.close() </div>
<br />
One can alternatively use the with formulation, which close the file automatically once used:<br />
<div style="background-color: #fce5cd;">
with open('datas.dat', 'r') as f2:</div>
<div style="background-color: #fce5cd;">
lines = f2.readlines()</div>
<i> </i><br />
<i>lines</i> is not usable as it is, as it contains the lines of the file in string form, without separating the 2 values...<br />
<br />
We now define <i>data</i> as a list, which will contain the lines of the file in the right format: <br />
<span style="background-color: #fce5cd;">data = []</span><br />
<br />
Now, we loop on the lines to copy the values, but after splitting the line: <br />
<span style="background-color: #fce5cd;">for line in lines:</span><br />
<span style="background-color: #fce5cd;"> p = line.split()</span><br />
<span style="background-color: #fce5cd;"> data.append(p)</span><br />
<br />
One can do the 4 previous commands in a single one, using list comprehension:<br />
<span style="background-color: #fce5cd;">data = [line.split() for line in lines]</span> <br />
<br />
Let's have a look at the first element of the list data:<br />
<div style="background-color: #fce5cd;">
print(data[0])</div>
<span style="background-color: #cfe2f3;">['0.000000', '0.000000']</span><br />
<br />
Ok, the value of <i>x</i> and <i>y</i> are not anymore in a single variable, but to really use them, we have to:<br />
<ol>
<li>transform the list into a numpy array, to allow operations on the values.</li>
<li>transform this into floating point.</li>
</ol>
The 2 operations can be performed one after the other:<br />
<br />
<span style="background-color: #fce5cd;">data2 = np.array(data)</span><br />
<span style="background-color: #fce5cd;">data2 = data2.astype(float)</span><br />
<br />
or in a single command (notice the f in <i>asfarray</i>, to transform into float):<br />
<span style="background-color: #fce5cd;">data2 = np.asfarray(data)</span><br />
<br />
The numpy arrays have a lot of properties and method. One tells us the shape of the array:<br />
<span style="background-color: #fce5cd;">print(data2.shape)</span><br />
<span style="background-color: #cfe2f3;">(100, 2)</span><br />
<br />
To access the distinct columns, use:<br />
<div style="background-color: #fce5cd;">
x = data2[:,0]</div>
<div style="background-color: #fce5cd;">
y = data2[:,1]</div>
<br />
As the array is a numpy object, one can perform operations on it:<br />
<div style="background-color: #fce5cd;">
z = (data2[:,0]**2 + data2[:,1]**2)**0.5</div>
<br />
Numpy includes functions to read ascii files and return float arrays directly. The following is then very compact and efficient:<br />
<span style="background-color: #fce5cd;">data3 = np.loadtxt('datas.dat')</span><br />
<span style="background-color: #fce5cd;"><br /></span><br />
An even more complete tool to read ascii file is the following:<br />
<span style="background-color: #fce5cd;">data4 = np.genfromtxt('datas.dat')</span><br />
It allows to read formated files, deal with undefined values, etc...<br />
One can even ask for this function to automatically put the result in 2 different variables:<br />
<div style="background-color: #fce5cd;">
x, y = np.loadtxt('datas.dat', unpack = True) </div>
<div style="background-color: #fce5cd;">
x, y = np.genfromtxt('datas.dat', unpack = True) </div>
<br />
Once you have the data2, data3, or data4, you can plot it (we'll see more on plots latter):<br />
<span style="background-color: #fce5cd;">plt.plot(data2[:,0], data2[:,1])</span><br />
<br />
A small comment on the order of the elements in arrays in Python:<br />
There is two ways arrays can be stored: row- or column major. It has a
direct impact on the way one has to loop on the arrays. IDL is like
Fortran (column major) and Python is like C (row major). It means that
in Python, as you move linearly through the memory of an array, the
second dimension (rightmost) changes the fastest, while in IDL the first
(leftmost) dimension changes the fastest.<br />
Consequence on the loop order in Python:<br />
<br />
<pre class="programlisting" style="background-color: #fce5cd;"><span class="py-statement">for</span> i <span class="py-statement">in</span> <span class="py-builtin">range</span>(100):
<span class="py-statement">for j in range(2):</span></pre>
<pre class="programlisting" style="background-color: #fce5cd;"><span class="py-statement"> ... data2[i,j] ...</span></pre>
<br />
Another comment on boundaries in arrays:<br />
<div style="background-color: #fce5cd;">
print(data2[0:4, 0])</div>
<span style="background-color: #cfe2f3;">[ 0. 0.0952 0.1904 0.285599]</span><br />
<i><span style="font-family: "Courier New",Courier,monospace;"></span></i>Some would expect the answer to be a 5 elements sub-array. The last index of 0:4 is not included. The rule is that a[n:m] return the elements between n and m-1.<br />
<br />
Something quite powerful: accessing the latest elements of an array:<br />
<div style="background-color: #fce5cd;">
print(data2[-1, :])</div>
<span style="background-color: #cfe2f3;">[ 9.424778 0. ] </span><br />
<div style="background-color: #fce5cd;">
print(data2[-5:-1, :])</div>
<div style="background-color: #cfe2f3;">
[[ 9.043979 0.371662]<br />
[ 9.139179 0.281733]<br />
[ 9.234378 0.189251]<br />
[ 9.329578 0.095056]]<br />
<i><span style="font-family: "Courier New",Courier,monospace;"></span></i></div>
<br />
The np.genfromtxt method allows to read formatted column tables, as for example the following file:<br />
<span style="font-size: x-small;"><span style="font-family: "Courier New",Courier,monospace;"># Line Iobs lambda relat_error Obs_code </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">H 1 4861A 1.00000 4861. 0.08000 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">H 1 6563A 2.8667 6563. 0.19467 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">H 1 4340A 0.4933 4340. 0.03307 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">H 1 4102A 0.2907 4102. 0.02229 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">H 1 3970A 0.1800 3970. 0.01253 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">N 2 6584A 2.1681 6584. 0.08686 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">N 2 121.7m 0.00446 1217000. 0.20000 Liu </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">O 1 6300A 0.0147 6300. 0.00325 Anabel </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">TOTL 2326A 0.07900 2326. 0.20000 Adams </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">C 2 157.6m 0.00856 1576000. 0.20000 Liu </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">O 1 63.17m 0.13647 631700. 0.10000 Liu </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">O 1 145.5m 0.00446 1455000. 0.200 Liu </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">TOTL 3727A 0.77609 3727. 0.200 Torres-Peimbert </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">S II 4070A 0.06174 4070. 0.200 Torres-Peimbert </span><br style="font-family: "Courier New",Courier,monospace;" />
<span style="font-family: "Courier New",Courier,monospace;">S II 4078A 0.06174 4078. 0.200 Torres-Peimbert </span><br style="font-family: "Courier New",Courier,monospace;" />
</span><br />
read using this command:<br />
<span style="background-color: #fce5cd;">obs = np.genfromtxt('observations.dat', skip_header=1,</span><br />
<span style="background-color: #fce5cd;"> dtype=["a11","float","float","float","a2","a2"], </span><br />
<span style="background-color: #fce5cd;"> delimiter=[11,8,9,8,2,2], </span><br />
<span style="background-color: #fce5cd;"> names = ['label', 'i_obs', 'lambda', 'e_obs', 'na', 'observer'],</span><br />
<span style="background-color: #fce5cd;"> usecols = (0,1,2,3,5)</span><br />
<span style="background-color: #fce5cd;"> )</span><br />
The skip_header tells not to read the first line. The dtype describes the type of the 6 columns. delimiter is the list of the sizes (in characters) of each column. The names are used to identify the columns. The usecols is used to specify which columns to output (here I don't need the 5th column, made of 2 spaces before the observer name)<br />
The data are accessible using for example:<br />
<span style="background-color: #fce5cd;">obs['i_obs']</span><br />
<br />
<i><span style="font-family: "Courier New",Courier,monospace;"></span></i><br />
To learn more in arrays indexing:<br />
<ul>
<li><a href="http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html">http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html</a> </li>
</ul>
To learn more about list comprehension:<br />
<ul>
<li><a href="http://docs.python.org/tutorial/datastructures.html#list-comprehensions">http://docs.python.org/tutorial/datastructures.html#list-comprehensions</a></li>
<li><a href="http://www.network-theory.co.uk/docs/pytut/ListComprehensions.html">http://www.network-theory.co.uk/docs/pytut/ListComprehensions.html </a></li>
<li><a href="http://www.secnetix.de/olli/Python/list_comprehensions.hawk">http://www.secnetix.de/olli/Python/list_comprehensions.hawk</a><br />
</li>
</ul>
More about <i>np.loadtxt, np.genfromtxt</i>:<br />
<ul>
<li><a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html">http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html</a></li>
<li><a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html">http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html</a><br />
</li>
</ul>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com0tag:blogger.com,1999:blog-8801264936720939312.post-79558982151868529842012-02-25T21:01:00.001-06:002014-09-03T16:23:50.783-05:00Read a 2 columns file and plot the result<div>
<h3>
<b>"Reading and writing files" is available in the 2014 version of the lecture here: </b><b><a href="http://python-astro.blogspot.mx/2014/09/interacting-with-files-reading-writing.html">http://python-astro.blogspot.mx/2014/09/interacting-with-files-reading-writing.html</a></b></h3>
<div style="color: blue;">
<b><br /></b></div>
<div style="color: blue;">
<b>In this second post, we will read the file we wrote in the previous one (datas.dat), and use the data to draw a plot.</b></div>
</div>
<br />
We will now have to import the plotting package, which is part of <i>matplotlib</i> (we will also need numpy):<br />
<br />
<span style="background-color: #fce5cd;">import matplotlib.pyplot as plt</span><br />
<span style="background-color: #fce5cd;">import numpy as np </span><br />
There is a lot of different ways to read a file, depending if it is ASCII or fits or Binary, if we want to extract only some columns, if we know the format of the data, etc... Even for a given data file, for example the one we did last post, I'll present various ways, to illustrate some functions of python.<br />
<br />
We now open the file with the 'r' parameter, to specify that we will read it:<br />
<span style="background-color: #fce5cd;">f2 = open('datas.dat', 'r') </span><br />
<i>'r'</i> is optional.<br />
<br />
We read the whole file in a single variable, we will deal with it latter. If the file is really big, and/or we only want a small part of it (we are looking for some line, or only one column is of interest) we may have to proceed differently:<br />
<div style="background-color: white;">
<span style="background-color: #fce5cd;">lines = f2.readlines()</span> </div>
<br />
Now that all the datas are in the <i>lines</i> variable, we can close the file and forget it:<br />
<span style="background-color: #fce5cd;">f2.close() </span><br />
<br />
<i>lines </i>contains the 200 values from the 2 vectors x and y, but we will have to extract this information, as it is now in an unusable form.<br />
First we can ask for the type of the variable:<br />
<div style="background-color: #fce5cd;">
type(lines)</div>
<span style="background-color: #cfe2f3;">list </span><br />
A list is... a list of objects. It can handle objects of different types, even other lists.<br />
<br />
We can print the length of the list, it is the number of elements it contains:<br />
<div style="background-color: #fce5cd;">
len(lines)</div>
<span style="background-color: #cfe2f3;">100</span><br />
<br />
We can print the first element:<br />
<div style="background-color: #fce5cd;">
print(lines[0])</div>
<span style="background-color: #cfe2f3;">0.000000 0.000000 </span> <br />
It seems it is the first line of our file, nice!<br />
<br />
But it is a single string...:<br />
<span style="background-color: #fce5cd;">print(type(lines[0]))</span><br />
<span style="background-color: #cfe2f3;"><type 'str'></span><br />
<br />
We will have to use a method of the object string to split the line into the 2 values:<br />
<span style="background-color: #fce5cd;">print(lines[0].split())</span><br />
<span style="background-color: #cfe2f3;">['0.000000', '0.000000']</span><br />
The result is another list, with 2 elements, each one being a string with the value of x and y corresponding to the line.<br />
<br />
We can now access to the value of x, and transform it into a floating point value:<br />
<span style="background-color: #fce5cd;">print(float(lines[0].split()[0]))</span><br />
<span style="background-color: #cfe2f3;">0.0</span><br />
<br />
We will now loop on all the lines of the file (they are stored in the <i>lines</i> variable) and do the same.<br />
We want to store the result in 2 lists. We first then need to define these lists, that will later be filled by the values of x and y.<br />
<div style="background-color: #fce5cd;">
x1 = []</div>
<span style="background-color: #fce5cd;">y1 = []</span><br />
<br />
The two empty lists are now ready. We can loop on files without to know the number of elements, python take care of this:<br />
<div style="background-color: #fce5cd;">
for line in lines:</div>
<span style="background-color: #fce5cd;"> p = line.split()</span><br />
<span style="background-color: #fce5cd;"> x1.append(float(p[0]))</span><br />
<span style="background-color: #fce5cd;"> y1.append(float(p[1]))</span><br />
<br />
We use <i>p </i>as an intermediate variable, to store the list of 2 strings (like ['0.000000', '0.000000'] for the first line).<br />
<i>append</i> is a method applied to the lists x1 and y1, to insert a new element to the end of the list.<br />
<br />
We now have the x and y vectors in the 2 variables x1 and y1.<br />
We can transform the two lists into <i>numpy</i> vectors, to have more flexibility on them: <br />
xv = np.array(x1)<br />
yv = np.array(y1)<br />
We can plot them:<br />
<div style="background-color: #fce5cd;">
plt.plot(xv, yv)</div>
<br />
To see the plot, one last command must be sent:<br />
<span style="background-color: #fce5cd;">plt.show()</span><br />
A window should pop-up with a nice sin(x) plot!<br />
<br />
All the previous commands can be stored in a file, called for example Ex2.py and execute by python.<br />
From the shell:<br />
<i>LINUX> python Ex2.py</i><br />
or from with python<br />
<div style="background-color: #fce5cd;">
import Ex2</div>
or from ipython<br />
<div style="background-color: #fce5cd;">
%run Ex2</div>
<br />
Here is the script of Ex2.py:<br />
<br />
<div style="background-color: #d9ead3;">
# Need to import the plotting package:<br />
import matplotlib.pyplot as plt<br />
<br />
# Read the file. <br />
f2 = open('datas.dat', 'r')<br />
# read the whole file into a single variable, which is a list of every row of the file.<br />
lines = f2.readlines()<br />
f2.close()<br />
<br />
# initialize some variable to be lists:<br />
x1 = []<br />
y1 = []<br />
<br />
# scan the rows of the file stored in lines, and put the values into some variables:<br />
for line in lines:<br />
p = line.split()<br />
x1.append(float(p[0]))<br />
y1.append(float(p[1]))<br />
<br />
xv = np.array(x1)<br />
yv = np.array(y1)<br />
<br />
<br />
# now, plot the data:<br />
plt.plot(xv, yv)<br />
<br />
plt.show()</div>
<div style="background-color: #d9ead3;">
<br /></div>
<br />
<br />
To learn more on lists:<br />
<ul>
<li><a href="http://docs.python.org/tutorial/introduction.html#lists">http://docs.python.org/tutorial/introduction.html#lists </a></li>
<li><a href="http://docs.python.org/tutorial/datastructures.html">http://docs.python.org/tutorial/datastructures.html</a></li>
<li><a href="http://www.tutorialspoint.com/python/python_lists.htm">http://www.tutorialspoint.com/python/python_lists.htm</a></li>
</ul>
To learn more on string manipulations:<br />
<ul>
<li><a href="http://docs.python.org/library/string.html">http://docs.python.org/library/string.html</a> </li>
</ul>
To learn more on plot (but we'll coma back to this latter):<br />
<ul>
<li><a href="http://matplotlib.sourceforge.net/users/pyplot_tutorial.html">http://matplotlib.sourceforge.net/users/pyplot_tutorial.html </a></li>
</ul>
Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com17tag:blogger.com,1999:blog-8801264936720939312.post-91337782113797475182012-02-25T15:54:00.001-06:002012-03-05T11:27:12.528-06:00Write vectors as an ASCII file<div style="color: blue;">
<b>In this first post we will create 2 vectors and write them into a file.</b></div>
<br />
Python can be used interactively, with classical call to python from the command line or by running iPython. In the following, I will mainly use iPython, which provides a lot of tools to help the user interacting with its variables and programs.<br />
<br />
Use copy-paste to execute the commands found in <span style="background-color: #fce5cd;">this color</span> in the following, in an interactive python or ipython shell. The outputs from python are <span style="background-color: #cfe2f3;">in this color.</span><br />
<br />
We can also put all the commands in a file and execute it from python (we'll come back to this latter).<br />
<br />
As we want to use vectors and perform operations on them, we need an additional library called <b>numpy.</b><br />
To import this library, we will write at the very beginning of the session (or the file):<br />
<span style="background-color: #fce5cd;">import numpy</span><i> </i><br />
<br />
Everything depending on numpy will be called usin<span style="background-color: white;">g</span><span style="background-color: white;"> numpy.</span><span style="background-color: white;"> before its name (Notice the point after numpy). For example, if we need to create an array with 10 elements named </span><span style="background-color: white;">a</span><span style="background-color: white;">, w</span>e do:<br />
<br />
<div style="background-color: #fce5cd;">
a = numpy.arange(10)</div>
<br />
Given that it's a little heavy to have to<span style="background-color: white;"> writ</span><span style="background-color: white;">e numpy.</span><span style="background-color: white;"> a l</span>ot of times, we can alias this when importing the package, so that we only need to use <b>np.</b>:<br />
<br />
<div style="background-color: #fce5cd;">
import numpy as np</div>
<div style="background-color: #fce5cd;">
a = np.arange(10)</div>
<br />
You can print out the value of <i>a</i> by just typing "a" and ENTER:<br />
<span style="background-color: #fce5cd;">a </span><br />
<span style="background-color: #cfe2f3;">array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) </span><br />
<br />
The output can differ if you use iPython:<br />
<span style="background-color: #fce5cd;">In [2]: a</span><br />
<span style="background-color: #cfe2f3;">Out[2]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])</span><br />
<br />
One can also print a:<br />
<div style="background-color: #fce5cd;">
print a</div>
<div style="background-color: #cfe2f3;">
[0 1 2 3 4 5 6 7 8 9]</div>
<br />
We want to create 2 vectors: <i>x</i> and <i>y = sin(x)</i>. We will first define some variables holding the parameters of the vectors:<br />
<div style="background-color: #fce5cd;">
x_min = 0.<br />
x_max = 3 * np.pi<br />
n_steps = 100</div>
<br />
To generate an array of n_steps values ranging from x_min to x_max, use linspace:<br />
<span style="background-color: #fce5cd;">x = np.linspace(x_min, x_max, n_steps)</span><br />
<br />
<i>x</i> is a numpy array, which support arithmetic operations like :<br />
<span style="background-color: #fce5cd;">y = x**2</span><br />
It can be used as argument for functions who will apply on each element of the array and return an array of the result. The trigonometrical functions are also available from numpy, so we use <i>np.</i> to access it:<br />
<span style="background-color: #fce5cd;">y = np.sin(x)</span><br />
<br />
We now want to store the 2 vectors as 2 columns of an ASCII file.<br />
We first open the file, named 'datas'. We tell python that we will write to it, using 'w':<br />
<span style="background-color: #fce5cd;">f = open('datas.dat', 'w')</span><br />
f is an object (everything in Python is an object ;-) ), and one can play with some function that apply to <i>f</i> or some variable related to <i>f</i> directly by calling them using <i>f.</i> (notice the point after f):<br />
<div style="background-color: #fce5cd;">
print f.name, f.mode</div>
<span style="background-color: #cfe2f3;">datas w</span><br />
We have to loop on the different elements of <i>x</i> and <i>y</i> to write them. The indexes in Python start at 0. The function <i>np.arange(n)</i> returns a list of integer from 0 to n-1, exactly what we need to scan the values of the two arrays <i>x</i> and <i>y</i>.<br />
The loop is classically defined by a <i>for</i> instruction. In Python, the start of a block (for, while, if, etc.. blobk) is delimited by :, and the block is defined by indentation:<br />
<span style="background-color: #fce5cd;">for i in np.arange(n_steps): </span><br />
<span style="background-color: #fce5cd;"> f.write('{0:f} {1:f} \n'.format(x[i], y[i]))</span><br />
<span style="background-color: #fce5cd;">f.close()</span><br />
<br />
We see that the f.write is indented: this is the block (here a single line, but can be very bigger) on which the loop is run. The instruction to close the file <i>f.close()</i> is written at the beginning of the line, no more indentation, so the loop is terminated.<br />
The format to write <i>x[i]</i> and <i>y[i]</i> is '{0:f} {1:f} \n'.<br />
{0:f} is to tell python to write the first variable in <i>format(x[i], y[i])</i>, namely x[i], as floating points.<br />
\n is to go to the next line. We will have a post dedicated to formats latter, you can already have a look at some links below.<br />
<br />
One can copy-paste all the previous commands in a python session, or write all in a file, to be run directly.<br />
<br />
The python program is there:<br />
<br />
<span style="background-color: #d9ead3;">#! /usr/bin/env python<br />#Need to import numpy:</span><br />
<span style="background-color: #d9ead3;">import numpy as np</span><br />
<br />
<span style="background-color: #d9ead3;"># define some variables:</span><br />
<span style="background-color: #d9ead3;">x_min = 0.</span><br />
<span style="background-color: #d9ead3;">x_max = 3 * np.pi</span><br />
<span style="background-color: #d9ead3;">n_steps = 100</span><br />
<span style="background-color: #d9ead3;"># create a vector from x_min to x_max, with n_steps elements</span><br />
<span style="background-color: #d9ead3;">x = np.linspace(x_min, x_max, n_steps)</span><br />
<span style="background-color: #d9ead3;"># create a 2nd vector, as a function of x</span><br />
<span style="background-color: #d9ead3;">y = np.sin(x)</span><br />
<br />
<span style="background-color: #d9ead3;"># Write x and y into a file</span><br />
<span style="background-color: #d9ead3;">f = open('datas.dat', 'w')</span><br />
<span style="background-color: #d9ead3;">for i in np.arange(n_steps):</span><br />
<span style="background-color: #d9ead3;"># could have been: for xx, yy in zip(x,y):</span><br />
<span style="background-color: #d9ead3;"> f.write('{0:f} {1:f} \n'.format(x[i], y[i])) </span><br />
<span style="background-color: #d9ead3;">f.close()</span><br />
<br />
Save it somewhere with the name Ex1.py.<br />
You can run it from the shell by:<br />
LINUX> python Ex1.py<br />
<br />
With the help of the first line <i>#! /usr/bin/env python</i>, once the Ex1.py file is set to executable (<i>chmod a+x Ex1.py</i>), it's possible to run it directly:<br />
LINUX> ./Ex1.py<br />
<br />
or from python:<br />
<span style="background-color: #fce5cd;">import Ex1</span><br />
Notice that import actually also executes any instruction in the imported file.<br />
<br />
or from ipython<br />
In [1]: %run Ex1<br />
<br />
The file "datas.dat" should contains the x and y values, in 2 columns, like this:<br />
0.000000 0.000000 <br />
0.095200 0.095056 <br />
0.190400 0.189251 <br />
0.285599 0.281733 <br />
0.380799 0.371662 <br />
0.475999 0.458227 <br />
0.571199 0.540641 <br />
0.666398 0.618159 <br />
0.761598 0.690079 <br />
0.856798 0.755750<br />
...<br />
<br />
To lear more on files:<br />
<ul>
<li><a href="http://docs.python.org/library/functions.html?highlight=open#open">http://docs.python.org/library/functions.html?highlight=open#open</a></li>
<li><a href="http://www.tutorialspoint.com/python/python_files_io.htm">http://www.tutorialspoint.com/python/python_files_io.htm</a></li>
</ul>
To learn more on string format:<br />
<ul>
<li>The % format: </li>
<ul>
<li><a href="http://docs.python.org/library/stdtypes.html#string-formatting-operations">http://docs.python.org/library/stdtypes.html#string-formatting-operations</a></li>
<li><a href="http://infohost.nmt.edu/tcc/help/pubs/lang/pytut/str-format.html">http://infohost.nmt.edu/tcc/help/pubs/lang/pytut/str-format.html</a></li>
</ul>
<li>The new format using {}: </li>
<ul>
<li><a href="http://www.python.org/dev/peps/pep-3101/">http://www.python.org/dev/peps/pep-3101/</a> </li>
<li><a href="http://knowledgestockpile.blogspot.com/2011/01/string-formatting-in-python_09.html">http://knowledgestockpile.blogspot.com/2011/01/string-formatting-in-python_09.html </a></li>
</ul>
</ul>
<ul>
<li>Comparing both: <a href="http://stackoverflow.com/questions/5082452/python-string-formatting-vs-format">http://stackoverflow.com/questions/5082452/python-string-formatting-vs-format </a></li>
</ul>Christophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com1tag:blogger.com,1999:blog-8801264936720939312.post-6278624933310598942012-02-25T11:36:00.000-06:002012-03-08T23:36:33.565-06:00BienvenidosHola!<br />
<br />
There is so many people who wants to follow the Python course at IA-UNAM that I decided to create a blog to share stuff. This way those who cannot assist the lecture can at least follow the advances and download the scripts and programs. I can also use this blog in real time while in the lecture.<br />
<br />
The useful packages that are recommended to install and some links are listed <a href="http://python-astro.blogspot.com/p/requiered-packages-and-important-links.html">THERE</a>. I will try
to keep this page uptodate, it is also accessible directly from the
right column of the blog.<br />
<br />
The messages are updated as we need it, so don't trust the publication date...<br />
<br />
So, bienvenidos, comments welcome ;-)<br />
ChristopheChristophehttp://www.blogger.com/profile/06332529829637156617noreply@blogger.com2