forked from python4astronomers/python4astronomers.github.com
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathnumpy_scipy.html
More file actions
975 lines (917 loc) · 65.4 KB
/
Copy pathnumpy_scipy.html
File metadata and controls
975 lines (917 loc) · 65.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>NumPy — Python4Astronomers 1.1 documentation</title>
<link rel="stylesheet" href="../_static/sphinxdoc.css" type="text/css" />
<link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: '../',
VERSION: '1.1',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="../_static/jquery.js"></script>
<script type="text/javascript" src="../_static/underscore.js"></script>
<script type="text/javascript" src="../_static/doctools.js"></script>
<link rel="top" title="Python4Astronomers 1.1 documentation" href="../index.html" />
<link rel="up" title="Core packages for analysis: IPython, NumPy, and SciPy" href="core.html" />
<link rel="next" title="CONTEST: Make a fun bouncing balls demo" href="../contest/bounce.html" />
<link rel="prev" title="IPython" href="ipython.html" />
<script type="text/javascript">
$(document).ready(function(){
$(".flip0").click(function(){
$(".panel0").slideToggle("normal");
});
$(".flip1").click(function(){
$(".panel1").slideToggle("normal");
});
$(".flip2").click(function(){
$(".panel2").slideToggle("normal");
});
$(".flip3").click(function(){
$(".panel3").slideToggle("normal");
});
$(".flip4").click(function(){
$(".panel4").slideToggle("normal");
});
$(".flip5").click(function(){
$(".panel5").slideToggle("normal");
});
$(".flip6").click(function(){
$(".panel6").slideToggle("normal");
});
$(".flip7").click(function(){
$(".panel7").slideToggle("normal");
});
$(".flip8").click(function(){
$(".panel8").slideToggle("normal");
});
$(".flip9").click(function(){
$(".panel9").slideToggle("normal");
});
});
</script>
<style type="text/css">
div.panel0,p.flip0
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip0
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel0
{
display:none;
}
div.panel1,p.flip1
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip1
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel1
{
display:none;
}
div.panel2,p.flip2
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip2
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel2
{
display:none;
}
div.panel3,p.flip3
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip3
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel3
{
display:none;
}
div.panel4,p.flip4
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip4
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel4
{
display:none;
}
div.panel5,p.flip5
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip5
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel5
{
display:none;
}
div.panel6,p.flip6
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip6
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel6
{
display:none;
}
div.panel7,p.flip7
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip7
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel7
{
display:none;
}
div.panel8,p.flip8
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip8
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel8
{
display:none;
}
div.panel9,p.flip9
{
font-size: 0.9em;
margin: 0;
padding: 0.1em 0 0.1em 0.5em;
border-bottom: 1px solid #86989B;
}
p.flip9
{
color: white;
font-weight: bold;
background-color: #AFC1C4;
}
div.panel9
{
display:none;
}
</style>
<script type="text/javascript">
var _gaq = _gaq || [];
_gaq.push(['_setAccount', 'UA-18709417-2']);
_gaq.push(['_trackPageview']);
(function() {
var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
})();
</script>
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="../contest/bounce.html" title="CONTEST: Make a fun bouncing balls demo"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="ipython.html" title="IPython"
accesskey="P">previous</a> |</li>
<li><a href="../index.html">Python4Astronomers 1.1 documentation</a> »</li>
<li><a href="core.html" accesskey="U">Core packages for analysis: IPython, NumPy, and SciPy</a> »</li>
</ul>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<h3><a href="../index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference internal" href="#">NumPy</a><ul>
<li><a class="reference internal" href="#setup">Setup</a></li>
<li><a class="reference internal" href="#read-in-the-2-d-image">Read in the 2-d image</a></li>
<li><a class="reference internal" href="#numpy-basics">NumPy basics</a></li>
<li><a class="reference internal" href="#plot-the-spatial-profile-and-raw-spectrum">Plot the spatial profile and raw spectrum</a></li>
<li><a class="reference internal" href="#filter-cosmic-rays-from-the-background">Filter cosmic rays from the background</a></li>
<li><a class="reference internal" href="#fit-the-background">Fit the background</a></li>
<li><a class="reference internal" href="#sum-the-source-signal">Sum the source signal</a></li>
</ul>
</li>
<li><a class="reference internal" href="#scipy">SciPy</a></li>
</ul>
<h4>Previous topic</h4>
<p class="topless"><a href="ipython.html"
title="previous chapter">IPython</a></p>
<h4>Next topic</h4>
<p class="topless"><a href="../contest/bounce.html"
title="next chapter">CONTEST: Make a fun bouncing balls demo</a></p>
<h3>This Page</h3>
<ul class="this-page-menu">
<li><a href="../_sources/core/numpy_scipy.txt"
rel="nofollow">Show Source</a></li>
</ul>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="../search.html" method="get">
<input type="text" name="q" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="numpy">
<h1>NumPy<a class="headerlink" href="#numpy" title="Permalink to this headline">¶</a></h1>
<p><a class="reference internal" href="#numpy">NumPy</a> is at the core of nearly every scientific Python application or
module since it provides a fast N-d array datatype that can be manipulated in a
vectorized form. This will be familiar to users of IDL or Matlab.</p>
<p>NumPy has a good and systematic <a class="reference external" href="http://www.scipy.org/Tentative_NumPy_Tutorial">basic tutorial</a> available. It is highly
recommended that you read this tutorial to fill in the gaps left by this
workshop, but on its own it’s a bit dry for the impatient astronomer.</p>
<p>Here we’ll learn NumPy by performing a very simple reduction of a
2-dimensional long slit spectrum (3C120 from HST/STIS):</p>
<ul class="simple">
<li>Read in the 2-d image</li>
<li>Plot the spatial profile and raw spectrum</li>
<li>Filter cosmic rays from the background</li>
<li>Fit for the background and subtract</li>
<li>Sum the source signal</li>
</ul>
<table border="1" class="docutils">
<colgroup>
<col width="51%" />
<col width="49%" />
</colgroup>
<thead valign="bottom">
<tr class="row-odd"><th class="head"><strong>2-d longslit image</strong></th>
<th class="head"><strong>Final 1-d spectrum</strong></th>
</tr>
</thead>
<tbody valign="top">
<tr class="row-even"><td><a class="first last reference internal image-reference" href="../_images/3c120.png"><img alt="../_images/3c120.png" src="../_images/3c120.png" style="width: 378.0px; height: 290.5px;" /></a>
</td>
<td><a class="first last reference internal image-reference" href="../_images/3c120_spec.gif"><img alt="../_images/3c120_spec.gif" src="../_images/3c120_spec.gif" style="width: 382.5px; height: 306.0px;" /></a>
</td>
</tr>
</tbody>
</table>
<div class="section" id="setup">
<h2>Setup<a class="headerlink" href="#setup" title="Permalink to this headline">¶</a></h2>
<p>Before going further you need to get the example data and script files for
the workshop. Now that you have a working Python installation we can do this
without worrying about details of the platform (e.g. linux has wget,
Mac has curl, Windows might not have tar, etc etc).</p>
<p>Now start IPython (“ipython –pylab”) or use your existing session and enter:</p>
<div class="highlight-python"><pre>import urllib2, tarfile
url = 'http://python4astronomers.github.com/core/core_examples.tar'
tarfile.open(fileobj=urllib2.urlopen(url), mode='r|').extractall()
cd py4ast/core
ls</pre>
</div>
<p>Leave this IPython session open for the rest of the workshop.</p>
<div class="admonition-exercise-for-the-interested-reader-how-did-that-code-above-work admonition">
<p class="first admonition-title">Exercise (for the interested reader): How did that code above work?</p>
<p class="last">Explain what’s happening in each part of the previous code snippet to grab
the file at a URL and untar it. Google on “python urllib2” and “python
tarfile” to find the relevant module docs. Figure out how you would
use the <tt class="docutils literal"><span class="pre">tarfile</span></tt> module to create a tarfile.</p>
</div>
<p class="flip0">Click to Show/Hide Solution</p> <div class="panel0"><ul class="simple">
<li><tt class="docutils literal"><span class="pre">urllib2.urlopen(url)</span></tt> opens the URL as a streaming file-like object</li>
<li><tt class="docutils literal"><span class="pre">mode='r|'</span> <span class="pre">means</span> <span class="pre">``tarfile</span></tt> is expecting a streaming file-like object
with no ability to seek in the file</li>
<li><tt class="docutils literal"><span class="pre">tarfile.open(..).extractall</span></tt> then extracts the tar archive</li>
</ul>
<p>Creating a tarfile is left for the reader to solve.</p>
</div></div>
<div class="section" id="read-in-the-2-d-image">
<h2>Read in the 2-d image<a class="headerlink" href="#read-in-the-2-d-image" title="Permalink to this headline">¶</a></h2>
<p>First read in the long-slit spectrum data. The standard file format available
for download from <a class="reference external" href="http://archive.stsci.edu/hst/">MAST</a> is a FITS file with
three identically sized images providing the 2-d spectral intensity, error
values, and data quality for each pixel. The slit direction is along the rows
(up and down) and wavelength is in columns (left to right).</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">astropy.io</span> <span class="kn">import</span> <span class="n">fits</span>
<span class="n">hdus</span> <span class="o">=</span> <span class="n">fits</span><span class="o">.</span><span class="n">open</span><span class="p">(</span><span class="s">'3c120_stis.fits.gz'</span><span class="p">)</span>
<span class="n">hdus</span>
</pre></div>
</div>
<p>Use the <tt class="docutils literal"><span class="pre">?</span></tt> to get a little more detail on the <tt class="docutils literal"><span class="pre">hdus</span></tt> object:</p>
<div class="highlight-python"><pre>hdus?</pre>
</div>
<p>Now give meaningful names to each of the three images that are available in the
FITS HDU list. You can access element <tt class="docutils literal"><span class="pre">n</span></tt> in a list with the index <tt class="docutils literal"><span class="pre">[n]</span></tt>,
where the count starts from 0:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">primary</span> <span class="o">=</span> <span class="n">hdus</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="c"># Primary (NULL) header data unit</span>
<span class="n">img</span> <span class="o">=</span> <span class="n">hdus</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="c"># Intensity data</span>
<span class="n">err</span> <span class="o">=</span> <span class="n">hdus</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="c"># Error per pixel</span>
<span class="n">dq</span> <span class="o">=</span> <span class="n">hdus</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="c"># Data quality per pixel</span>
</pre></div>
</div>
<p>Next have a look at the images using one of the standard Matplotlib plotting
functions:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">img</span><span class="p">)</span>
</pre></div>
</div>
<p>As you can see, it is hard to see things. So, let’s set a few option for this
plot. First, we want the origin in the lower left instead of the upper left
corner:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">img</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">)</span>
</pre></div>
</div>
<p>Second, let’s change the scaling to something more sensible. By default,
<tt class="docutils literal"><span class="pre">plt.imshow()</span></tt> scales the colorbar from the minimum to the maximum value. In
our case that is not the best option. We can set a lower and upper bound and
add a colorbar to our plot:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">img</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">,</span> <span class="n">vmin</span> <span class="o">=</span> <span class="o">-</span><span class="mi">10</span><span class="p">,</span> <span class="n">vmax</span> <span class="o">=</span> <span class="mi">65</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">colorbar</span><span class="p">()</span>
</pre></div>
</div>
<p>Your plot should not look like this (it is possible that the colormap differs,
if your matplotlib has different defaults set).</p>
<a class="reference internal image-reference" href="../_images/imgview_img.png"><img alt="../_images/imgview_img.png" src="../_images/imgview_img.png" style="width: 406.0px; height: 306.0px;" /></a>
<div class="admonition-exercise-view-the-error-and-data-quality-images admonition">
<p class="first admonition-title">Exercise: View the error and data quality images</p>
<p class="last">Bring up a viewer window for the other two images. Play with the toolbar
buttons on the lower-left (hint: try the four on the right first, then
imagine a web browser for the three on the left). Does the save button
work for you?</p>
</div>
<p class="flip1">Click to Show/Hide Solution</p> <div class="panel1"><div class="highlight-python"><div class="highlight"><pre><span class="c"># Errors</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">err</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">,</span> <span class="n">vmin</span> <span class="o">=</span> <span class="mi">5</span><span class="p">,</span> <span class="n">vmax</span> <span class="o">=</span> <span class="mi">25</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">colorbar</span><span class="p">()</span>
<span class="c"># Data quality</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">dq</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">,</span> <span class="n">vmax</span> <span class="o">=</span> <span class="mi">25</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">colorbar</span><span class="p">()</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/imgview_err.png"><img alt="../_images/imgview_err.png" src="../_images/imgview_err.png" style="width: 406.0px; height: 306.0px;" /></a>
<a class="reference internal image-reference" href="../_images/imgview_dq.png"><img alt="../_images/imgview_dq.png" src="../_images/imgview_dq.png" style="width: 406.0px; height: 306.0px;" /></a>
</div><p>Now discover a little bit about the images you have read in, first with <tt class="docutils literal"><span class="pre">?</span></tt>:</p>
<div class="highlight-python"><pre>img?</pre>
</div>
<p>Next use <tt class="docutils literal"><span class="pre">help</span></tt> and note the slightly different information that you get:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">help</span><span class="p">(</span><span class="n">img</span><span class="p">)</span>
</pre></div>
</div>
<p>Use tab completion to see all the methods in short form:</p>
<div class="highlight-python"><pre>img.<TAB></pre>
</div>
<p>Finally find the shape of the image and its minimum value:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">img</span><span class="o">.</span><span class="n">shape</span> <span class="c"># Get the shape of img</span>
<span class="n">img</span><span class="o">.</span><span class="n">min</span><span class="p">()</span> <span class="c"># Call object method min with no arguments</span>
</pre></div>
</div>
</div>
<div class="section" id="numpy-basics">
<h2>NumPy basics<a class="headerlink" href="#numpy-basics" title="Permalink to this headline">¶</a></h2>
<p>Before going further on the spectral extraction project we need to learn about
a few key features of NumPy.</p>
<div class="section" id="making-arrays">
<h3>Making arrays<a class="headerlink" href="#making-arrays" title="Permalink to this headline">¶</a></h3>
<p>Arrays can be created in different ways. The “>>>” indicates the input to Python:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">10</span><span class="p">,</span> <span class="mi">20</span><span class="p">,</span> <span class="mi">30</span><span class="p">,</span> <span class="mi">40</span><span class="p">])</span> <span class="c"># create an array from a list of values</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([10, 20, 30, 40]</span>
<span class="gp">>>> </span><span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span> <span class="c"># create an array of 4 integers, from 0 to 3</span>
<span class="gp">>>> </span><span class="n">b</span>
<span class="go">array([0, 1, 2, 3]),</span>
<span class="gp">>>> </span><span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">10.0</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">)</span> <span class="c"># create a float array from 0 to 100 stepping by 0.1</span>
<span class="go">array([ 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ,</span>
<span class="go"> 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1,</span>
<span class="go"> 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2,</span>
<span class="go"> 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3,</span>
<span class="go"> 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5. , 5.1, 5.2, 5.3, 5.4,</span>
<span class="go"> 5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1, 6.2, 6.3, 6.4, 6.5,</span>
<span class="go"> 6.6, 6.7, 6.8, 6.9, 7. , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6,</span>
<span class="go"> 7.7, 7.8, 7.9, 8. , 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,</span>
<span class="go"> 8.8, 8.9, 9. , 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8,</span>
<span class="go"> 9.9]),</span>
<span class="gp">>>> </span><span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="mi">5</span><span class="p">)</span> <span class="c"># create an array of 5 evenly spaced samples from -pi to pi</span>
<span class="go">array([-3.14159265, -1.57079633, 0. , 1.57079633, 3.14159265]))</span>
</pre></div>
</div>
<p>New arrays can be obtained by operating with existing arrays:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">+</span> <span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="c"># elementwise operations</span>
<span class="go">array([10, 21, 34, 49])</span>
</pre></div>
</div>
<p>Arrays may have more than one dimension:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">f</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">([</span><span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">])</span> <span class="c"># 3 x 4 float array of ones</span>
<span class="gp">>>> </span><span class="n">f</span>
<span class="go">array([[ 1., 1., 1., 1.],</span>
<span class="go"> [ 1., 1., 1., 1.],</span>
<span class="go"> [ 1., 1., 1., 1.]]),</span>
<span class="gp">>>> </span><span class="n">g</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">],</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">int</span><span class="p">)</span> <span class="c"># 2 x 3 x 4 int array of zeros</span>
<span class="go">array([[[0, 0, 0, 0],</span>
<span class="go"> [0, 0, 0, 0],</span>
<span class="go"> [0, 0, 0, 0]],</span>
<span class="go"> [[0, 0, 0, 0],</span>
<span class="go"> [0, 0, 0, 0],</span>
<span class="go"> [0, 0, 0, 0]]]),</span>
<span class="gp">>>> </span><span class="n">i</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros_like</span><span class="p">(</span><span class="n">f</span><span class="p">)</span> <span class="c"># array of zeros with same shape/type as f</span>
<span class="go">array([[ 0., 0., 0., 0.],</span>
<span class="go"> [ 0., 0., 0., 0.],</span>
<span class="go"> [ 0., 0., 0., 0.]]))</span>
</pre></div>
</div>
<p>You can change the dimensions of existing arrays:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">w</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">12</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">w</span><span class="o">.</span><span class="n">shape</span> <span class="o">=</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span> <span class="mi">4</span><span class="p">]</span> <span class="c"># does not modify the total number of elements</span>
<span class="go">array([[ 0, 1, 2, 3],</span>
<span class="go"> [ 4, 5, 6, 7],</span>
<span class="go"> [ 8, 9, 10, 11]]),</span>
<span class="gp">>>> </span><span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">5</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">x</span>
<span class="go">array([0, 1, 2, 3, 4]),</span>
<span class="gp">>>> </span><span class="n">y</span> <span class="o">=</span> <span class="n">x</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">5</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">y</span> <span class="o">=</span> <span class="n">x</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span> <span class="c"># Same thing but NumPy figures out correct length</span>
<span class="gp">>>> </span><span class="n">y</span>
<span class="go">array([[0],</span>
<span class="go"> [1],</span>
<span class="go"> [2],</span>
<span class="go"> [3],</span>
<span class="go"> [4]]))</span>
</pre></div>
</div>
<p>It is possible to operate with arrays of different dimensions as long
as they fit well (this is known as
<a class="reference external" href="http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html">broadcasting</a>
in NumPy):</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">x</span> <span class="o">+</span> <span class="n">y</span> <span class="o">*</span> <span class="mi">10</span>
<span class="go">array([[ 0, 1, 2, 3, 4],</span>
<span class="go"> [10, 11, 12, 13, 14],</span>
<span class="go"> [20, 21, 22, 23, 24],</span>
<span class="go"> [30, 31, 32, 33, 34],</span>
<span class="go"> [40, 41, 42, 43, 44]])</span>
</pre></div>
</div>
<div class="admonition-exercise-make-a-ripple admonition">
<p class="first admonition-title">Exercise: Make a ripple</p>
<p class="last">Calculate a surface <tt class="docutils literal"><span class="pre">z</span> <span class="pre">=</span> <span class="pre">cos(r)</span> <span class="pre">/</span> <span class="pre">(r</span> <span class="pre">+</span> <span class="pre">5)</span></tt> where <tt class="docutils literal"><span class="pre">r</span> <span class="pre">=</span> <span class="pre">sqrt(x**2</span> <span class="pre">+</span>
<span class="pre">y**2)</span></tt>. Set <tt class="docutils literal"><span class="pre">x</span></tt> to an array that goes from -20 to 20 stepping by 0.25
Make <tt class="docutils literal"><span class="pre">y</span></tt> the same as <tt class="docutils literal"><span class="pre">x</span></tt> but “transposed” using the <tt class="docutils literal"><span class="pre">reshape</span></tt> trick above.
Use <cite>plt.imshow</cite> to display the image of <tt class="docutils literal"><span class="pre">z</span></tt>.</p>
</div>
<p class="flip3">Click to Show/Hide Solution</p> <div class="panel3"><div class="highlight-python"><pre>x = np.arange(-20, 20, 0.25)
y = x.reshape(-1, 1)
r = np.sqrt(x**2 + y**2)
z = np.cos(r) / (r + 5)
plt.imshow(z, origin = 'lower)</pre>
</div>
<a class="reference internal image-reference" href="../_images/ripple.png"><img alt="../_images/ripple.png" src="../_images/ripple.png" style="width: 406.0px; height: 306.0px;" /></a>
</div></div>
<div class="section" id="array-access-and-slicing">
<h3>Array access and slicing<a class="headerlink" href="#array-access-and-slicing" title="Permalink to this headline">¶</a></h3>
<p>NumPy provides powerful methods for accessing array elements or particular subsets of an array,
e.g. the 4th column or every other row. This is called slicing. The outputs
below illustrate basic slicing, but you don’t need to type these examples:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">20</span><span class="p">)</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">4</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([[ 0, 1, 2, 3, 4],</span>
<span class="go"> [ 5, 6, 7, 8, 9],</span>
<span class="go"> [10, 11, 12, 13, 14],</span>
<span class="go"> [15, 16, 17, 18, 19]])</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">]</span> <span class="c"># select element in row 2, col 3 (counting from 0)</span>
<span class="go">13</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="p">:]</span> <span class="c"># select every element in row 2</span>
<span class="go">array([10, 11, 12, 13, 14])</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span> <span class="c"># select every element in col 0</span>
<span class="go">array([ 0, 5, 10, 15])</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="mi">0</span><span class="p">:</span><span class="mi">3</span><span class="p">,</span> <span class="mi">1</span><span class="p">:</span><span class="mi">3</span><span class="p">]</span>
<span class="go">array([[ 1, 2],</span>
<span class="go"> [ 6, 7],</span>
<span class="go"> [11, 12]])</span>
</pre></div>
</div>
<p>As a first practical
example plot column 300 of the longslit image to look at the spatial profile:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span> <span class="c"># Clear the existing plot -- by default matplotlib overplots.</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">img</span><span class="p">[:,</span> <span class="mi">300</span><span class="p">])</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/img_col300.png"><img alt="../_images/img_col300.png" src="../_images/img_col300.png" style="width: 406.0px; height: 306.0px;" /></a>
<p>The full slicing syntax also allows for a step size:</p>
<div class="highlight-python"><pre><slice> = i0:i1:step
array[<slice0>, <slice1>, ...]</pre>
</div>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">i0</span></tt> is the first index value (default is zero if not provided)</li>
<li><tt class="docutils literal"><span class="pre">i1</span></tt> is the index upper bound (default is last element index + 1)</li>
<li><tt class="docutils literal"><span class="pre">step</span></tt> is the step size (default is one). When <tt class="docutils literal"><span class="pre">step</span></tt> is not specified then the final ”:” is not required.</li>
</ul>
<div class="admonition-exercise-slice-the-error-array admonition">
<p class="first admonition-title">Exercise: Slice the error array</p>
<ul class="last simple">
<li>For row 254 of the error array <tt class="docutils literal"><span class="pre">err</span></tt> plot columns 10 to 200 stepping by 3.</li>
<li>Print a rectangular region slice of the data quality with rows 251 to 253 (inclusive) and columns 101 to
104 (inclusive). What did you learn about the index upper bound value?</li>
</ul>
</div>
<p class="flip2">Click to Show/Hide Solution</p> <div class="panel2"><div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">err</span><span class="p">[</span><span class="mi">254</span><span class="p">,</span> <span class="mi">10</span><span class="p">:</span><span class="mi">200</span><span class="p">:</span><span class="mi">3</span><span class="p">])</span>
<span class="n">dq</span><span class="p">[</span><span class="mi">251</span><span class="p">:</span><span class="mi">254</span><span class="p">,</span> <span class="mi">101</span><span class="p">:</span><span class="mi">105</span><span class="p">]</span>
</pre></div>
</div>
<p>The index upper bound <tt class="docutils literal"><span class="pre">i1</span></tt> is one more than the final index that gets
included in the slice. In other words the slice includes everything up to,
<em>but not including</em>, the index upper bound <tt class="docutils literal"><span class="pre">i1</span></tt>. There are good reasons for
this, but for now just accept and learn it.</p>
<a class="reference internal image-reference" href="../_images/err_row254.png"><img alt="../_images/err_row254.png" src="../_images/err_row254.png" style="width: 406.0px; height: 306.0px;" /></a>
</div></div>
</div>
<div class="section" id="plot-the-spatial-profile-and-raw-spectrum">
<h2>Plot the spatial profile and raw spectrum<a class="headerlink" href="#plot-the-spatial-profile-and-raw-spectrum" title="Permalink to this headline">¶</a></h2>
<p>Plot the spatial profile by summing along the wavelength direction:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">profile</span> <span class="o">=</span> <span class="n">img</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">profile</span><span class="p">)</span>
</pre></div>
</div>
<p>Now plot the spectrum by summing along the spatial direction:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">spectrum</span> <span class="o">=</span> <span class="n">img</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">spectrum</span><span class="p">)</span>
</pre></div>
</div>
<p>Since most of the sum is in the background region there is a lot of noise and
cosmic-ray contamination.</p>
<a class="reference internal image-reference" href="../_images/profile.png"><img alt="../_images/profile.png" src="../_images/profile.png" style="width: 406.0px; height: 306.0px;" /></a>
<a class="reference internal image-reference" href="../_images/spectrum_noisy.png"><img alt="../_images/spectrum_noisy.png" src="../_images/spectrum_noisy.png" style="width: 406.0px; height: 306.0px;" /></a>
<div class="admonition-exercise-use-slicing-to-make-a-better-spectrum-plot admonition">
<p class="first admonition-title">Exercise: Use slicing to make a better spectrum plot</p>
<p class="last">Use slicing to do the spectrum sum using only the rows in the image where
there is a signal from the source.
Hint: zoom into the profile plot to find the right row range.</p>
</div>
<p class="flip4">Click to Show/Hide Solution</p> <div class="panel4"><div class="highlight-python"><div class="highlight"><pre><span class="n">spectrum</span> <span class="o">=</span> <span class="n">img</span><span class="p">[</span><span class="mi">250</span><span class="p">:</span><span class="mi">260</span><span class="p">,</span> <span class="p">:]</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">spectrum</span><span class="p">)</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/spectrum_clean.png"><img alt="../_images/spectrum_clean.png" src="../_images/spectrum_clean.png" style="width: 406.0px; height: 306.0px;" /></a>
</div></div>
<div class="section" id="filter-cosmic-rays-from-the-background">
<h2>Filter cosmic rays from the background<a class="headerlink" href="#filter-cosmic-rays-from-the-background" title="Permalink to this headline">¶</a></h2>
<p>Plot five columns (wavelength) from the spectrum image as follows:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">img</span><span class="p">[:,</span> <span class="mi">254</span><span class="p">:</span><span class="mi">259</span><span class="p">])</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/img_row254_noisy.png"><img alt="../_images/img_row254_noisy.png" src="../_images/img_row254_noisy.png" style="width: 406.0px; height: 306.0px;" /></a>
<p>The basic idea in spectral extraction is to subtract out the background and sum
over rows with the source signal.</p>
<p>It’s evident that there are significant cosmic ray defects in the data. In
order to do a good job of subtracting the background we need to filter them
out. Doing this correctly in general is difficult and in reality one would
just use the answers already provided by STSci.</p>
<p><strong>Strategy</strong>: Use a median filter to smooth out single-pixel deviations. Then
use sigma-clipping to remove large variations between the actual and smoothed
image.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">scipy.signal</span>
<span class="n">img_sm</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">signal</span><span class="o">.</span><span class="n">medfilt</span><span class="p">(</span><span class="n">img</span><span class="p">,</span> <span class="mi">5</span><span class="p">)</span>
<span class="n">sigma</span> <span class="o">=</span> <span class="n">median</span><span class="p">(</span><span class="n">err</span><span class="p">)</span>
<span class="n">bad</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">img</span> <span class="o">-</span> <span class="n">img_sm</span><span class="p">)</span> <span class="o">/</span> <span class="n">sigma</span> <span class="o">></span> <span class="mf">8.0</span>
<span class="n">img_cr</span> <span class="o">=</span> <span class="n">img</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
<span class="n">img_cr</span><span class="p">[</span><span class="n">bad</span><span class="p">]</span> <span class="o">=</span> <span class="n">img_sm</span><span class="p">[</span><span class="n">bad</span><span class="p">]</span>
<span class="n">img_cr</span><span class="p">[</span><span class="mi">230</span><span class="p">:</span><span class="mi">280</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">img</span><span class="p">[</span><span class="mi">230</span><span class="p">:</span><span class="mi">280</span><span class="p">,:]</span> <span class="c"># Filter only for background</span>
</pre></div>
</div>
<p>Check if it worked:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">img_cr</span><span class="p">[:,</span> <span class="mi">254</span><span class="p">:</span><span class="mi">259</span><span class="p">])</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/img_row254_clean.png"><img alt="../_images/img_row254_clean.png" src="../_images/img_row254_clean.png" style="width: 406.0px; height: 306.0px;" /></a>
<p>This introduces the important concept of slicing with a <strong>boolean mask</strong>. Let’s
look at a smaller example:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">1</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">5</span><span class="p">])</span>
<span class="gp">>>> </span><span class="n">neg</span> <span class="o">=</span> <span class="p">(</span><span class="n">a</span> <span class="o"><</span> <span class="mi">0</span><span class="p">)</span> <span class="c"># Parentheses here for clarity but are not required</span>
<span class="gp">>>> </span><span class="n">neg</span>
<span class="go">array([False, False, True, False, True], dtype=bool)</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="n">neg</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([1, 4, 0, 4, 0])</span>
</pre></div>
</div>
<p>A slightly more complex example shows that this works the same on N-d arrays
and that you can compose logical expressions:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">25</span><span class="p">)</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">ok</span> <span class="o">=</span> <span class="p">(</span><span class="n">a</span> <span class="o">></span> <span class="mi">6</span><span class="p">)</span> <span class="o">&</span> <span class="p">(</span><span class="n">a</span> <span class="o"><</span> <span class="mi">17</span><span class="p">)</span> <span class="c"># "ok = a > 6 & a < 17" will FAIL!</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="o">~</span><span class="n">ok</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span> <span class="c"># Note the "logical not" operator</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([[ 0, 0, 0, 0, 0],</span>
<span class="go"> [ 0, 0, 7, 8, 9],</span>
<span class="go"> [10, 11, 12, 13, 14],</span>
<span class="go"> [15, 16, 0, 0, 0],</span>
<span class="go"> [ 0, 0, 0, 0, 0]])</span>
</pre></div>
</div>
<div class="admonition-exercise-intermediate-circular-region-slicing admonition">
<p class="first admonition-title">Exercise [intermediate]: circular region slicing</p>
<p class="last">Remember the surface <tt class="docutils literal"><span class="pre">z</span> <span class="pre">=</span> <span class="pre">cos(r)</span> <span class="pre">/</span> <span class="pre">(r</span> <span class="pre">+</span> <span class="pre">5)</span></tt> that you made previously. Set
<tt class="docutils literal"><span class="pre">z</span> <span class="pre">=</span> <span class="pre">0</span></tt> for every pixel of <tt class="docutils literal"><span class="pre">z</span></tt> that is within 10 units of (x,y) = (10, 15).</p>
</div>
<p class="flip5">Click to Show/Hide Solution</p> <div class="panel5"><div class="highlight-python"><div class="highlight"><pre><span class="n">dist</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">x</span><span class="o">-</span><span class="mi">10</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="p">(</span><span class="n">y</span><span class="o">-</span><span class="mi">15</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
<span class="n">mask</span> <span class="o">=</span> <span class="n">dist</span> <span class="o"><</span> <span class="mi">10</span>
<span class="n">z</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">z</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">)</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/ripple_masked.png"><img alt="../_images/ripple_masked.png" src="../_images/ripple_masked.png" style="width: 406.0px; height: 306.0px;" /></a>
</div><div class="admonition-detour-copy-versus-reference admonition">
<p class="first admonition-title">Detour: copy versus reference</p>
<dl class="docutils">
<dt><strong>Question</strong></dt>
<dd>In the median filtering commands above we wrote <tt class="docutils literal"><span class="pre">img_cr</span> <span class="pre">=</span> <span class="pre">img.copy()</span></tt>. Why
was that needed instead of just <tt class="docutils literal"><span class="pre">img_cr</span> <span class="pre">=</span> <span class="pre">img</span></tt>?</dd>
<dt><strong>Answer</strong></dt>
<dd>Because the statement <tt class="docutils literal"><span class="pre">img_cr</span> <span class="pre">=</span> <span class="pre">img</span></tt> would just create another reference
pointing to the underlying N-d array object that <tt class="docutils literal"><span class="pre">img</span></tt> references.</dd>
</dl>
<p>Variable names in Python are just pointers to the actual Python
object. To see this clearly do the following:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">8</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">b</span> <span class="o">=</span> <span class="n">a</span>
<span class="gp">>>> </span><span class="nb">id</span><span class="p">(</span><span class="n">a</span><span class="p">)</span> <span class="c"># Unique identifier for the object referred to by "a": arange(8)</span>
<span class="go">122333200</span>
<span class="gp">>>> </span><span class="nb">id</span><span class="p">(</span><span class="n">b</span><span class="p">)</span> <span class="c"># Unique identifier for the object referred to by "b": same ^^</span>
<span class="go">122333200</span>
<span class="gp">>>> </span><span class="n">b</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="mi">10</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([ 0, 1, 2, -10, 4, 5, 6, 7])</span>
</pre></div>
</div>
<p>After getting over the initial confusion this behavior is actually a good
thing because it is efficient and consistent within Python. If you really
need a copy of an array then use the copy() method as shown.</p>
<p><strong>BEWARE</strong> of one common pitfall: NumPy “basic” slicing like <tt class="docutils literal"><span class="pre">a[3:6]</span></tt>
does NOT make a copy:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">b</span> <span class="o">=</span> <span class="n">a</span><span class="p">[</span><span class="mi">3</span><span class="p">:</span><span class="mi">6</span><span class="p">]</span>
<span class="gp">>>> </span><span class="n">b</span>
<span class="go">array([-10, 4, 5])</span>
<span class="gp">>>> </span><span class="n">b</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mi">100</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([ 0, 1, 2, -10, 100, 5, 6, 7])</span>
</pre></div>
</div>
<p>However if you do arithmetic or boolean mask then a copy is always made:</p>
<div class="last highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
<span class="gp">>>> </span><span class="n">b</span> <span class="o">=</span> <span class="n">a</span><span class="o">**</span><span class="mi">2</span>
<span class="gp">>>> </span><span class="n">a</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mi">100</span>
<span class="gp">>>> </span><span class="n">a</span>
<span class="go">array([ 0, 100, 2, 3])</span>
<span class="gp">>>> </span><span class="n">b</span> <span class="c"># Still as expected after changing "a"</span>
<span class="go">array([0, 1, 4, 9])</span>
</pre></div>
</div>
</div>
</div>
<div class="section" id="fit-the-background">
<h2>Fit the background<a class="headerlink" href="#fit-the-background" title="Permalink to this headline">¶</a></h2>
<p>To subtract the background signal from the source region we want to fit a
quadratic to the background pixels and subtract that quadratic from the entire
image which includes the source region.</p>
<p>Let’s tackle a simpler problem first and fit the background for a single column:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">x</span> <span class="o">=</span> <span class="n">append</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">10</span><span class="p">,</span> <span class="mi">200</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mi">300</span><span class="p">,</span> <span class="mi">480</span><span class="p">))</span> <span class="c"># Background rows</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">img_cr</span><span class="p">[</span><span class="n">x</span><span class="p">,</span> <span class="mi">10</span><span class="p">]</span> <span class="c"># Background rows of column 10 of cleaned image</span>
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span>
<span class="n">pfit</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">polyfit</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="mi">2</span><span class="p">)</span> <span class="c"># Fit a 2nd order polynomial to (x, y) data</span>
<span class="n">yfit</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">polyval</span><span class="p">(</span><span class="n">pfit</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span> <span class="c"># Evaluate the polynomial at x</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">yfit</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">()</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/bkg_fit0.png"><img alt="../_images/bkg_fit0.png" src="../_images/bkg_fit0.png" style="width: 406.0px; height: 306.0px;" /></a>
<p>Now do this for every column and store the results in a background image:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">xrows</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="n">img_cr</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="c"># Array from 0 .. N_rows-1</span>
<span class="n">bkg</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros_like</span><span class="p">(</span><span class="n">img_cr</span><span class="p">)</span> <span class="c"># Empty image for background fits</span>
<span class="k">for</span> <span class="n">col</span> <span class="ow">in</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="n">img_cr</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]):</span> <span class="c"># Iterate over columns</span>
<span class="n">pfit</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">polyfit</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">img_cr</span><span class="p">[</span><span class="n">x</span><span class="p">,</span> <span class="n">col</span><span class="p">],</span> <span class="mi">2</span><span class="p">)</span> <span class="c"># Fit poly over bkg rows for col</span>
<span class="n">bkg</span><span class="p">[:,</span> <span class="n">col</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">polyval</span><span class="p">(</span><span class="n">pfit</span><span class="p">,</span> <span class="n">xrows</span><span class="p">)</span> <span class="c"># Eval poly at ALL row positions</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">bkg</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">,</span> <span class="n">vmin</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">vmax</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">colorbar</span><span class="p">()</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/bkg_fit1.png"><img alt="../_images/bkg_fit1.png" src="../_images/bkg_fit1.png" style="width: 406.0px; height: 306.0px;" /></a>
<p>Finally subtract this background and see if it worked:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">img_bkg</span> <span class="o">=</span> <span class="n">img_cr</span> <span class="o">-</span> <span class="n">bkg</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">img_bkg</span><span class="p">,</span> <span class="n">origin</span> <span class="o">=</span> <span class="s">'lower'</span><span class="p">,</span> <span class="n">vmin</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">vmax</span><span class="o">=</span><span class="mi">60</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">colorbar</span><span class="p">()</span>
</pre></div>
</div>
<table border="1" class="docutils">
<colgroup>
<col width="51%" />
<col width="49%" />
</colgroup>
<thead valign="bottom">
<tr class="row-odd"><th class="head"><strong>Background subtracted</strong></th>
<th class="head"><strong>Original</strong></th>
</tr>
</thead>
<tbody valign="top">
<tr class="row-even"><td><a class="first last reference internal image-reference" href="../_images/bkg_fit2.png"><img alt="../_images/bkg_fit2.png" src="../_images/bkg_fit2.png" style="width: 406.0px; height: 306.0px;" /></a>
</td>
<td><a class="first last reference internal image-reference" href="../_images/imgview_img.png"><img alt="../_images/imgview_img.png" src="../_images/imgview_img.png" style="width: 406.0px; height: 306.0px;" /></a>
</td>
</tr>
</tbody>
</table>
<div class="admonition-detour-vector-operations-versus-looping admonition">
<p class="first admonition-title">Detour: vector operations versus looping</p>
<p>If you are used to C or Fortran you might be wondering why jump through these
hoops with slicing and making sure everything is vectorized. The answer is
that pure Python is an interpreted dynamic language and hence doing loops is
<em>slow</em>. Try the following:</p>
<div class="highlight-python"><pre>size = 500000
x = np.arange(size)
a = np.zeros(size)
time for i in x: a[i] = x[i] / 2.0</pre>
</div>
<p>Now compare to the vectorized NumPy solution:</p>
<div class="highlight-python"><pre>x = np.arange(size)
time a = x / 2</pre>
</div>
<p class="last">Sometimes doing things in a vectorized way is not possible or just too
confusing. There is an art here and the basic answer is that if it runs
fast enough then you are good to go. Otherwise things need to be vectorized
or maybe coded in C or Fortran.</p>
</div>
</div>
<div class="section" id="sum-the-source-signal">
<h2>Sum the source signal<a class="headerlink" href="#sum-the-source-signal" title="Permalink to this headline">¶</a></h2>
<p>Now the final step is easy and is left as an exercise.</p>
<table border="1" class="docutils">
<colgroup>
<col width="51%" />
<col width="49%" />
</colgroup>
<thead valign="bottom">
<tr class="row-odd"><th class="head"><strong>Python for Astronomers Spectrum</strong></th>
<th class="head"><strong>HST official spectrum</strong></th>
</tr>
</thead>
<tbody valign="top">
<tr class="row-even"><td><a class="first last reference internal image-reference" href="../_images/spectrum_final.png"><img alt="../_images/spectrum_final.png" src="../_images/spectrum_final.png" style="width: 406.0px; height: 306.0px;" /></a>
</td>
<td><a class="first last reference internal image-reference" href="../_images/3c120_spec.gif"><img alt="../_images/3c120_spec.gif" src="../_images/3c120_spec.gif" style="width: 382.5px; height: 306.0px;" /></a>
</td>
</tr>
</tbody>
</table>
<div class="admonition-exercise-make-the-final-spectrum admonition">
<p class="first admonition-title">Exercise: Make the final spectrum</p>
<p class="last">Sum the rows of the background subtracted spectrum and plot. Hint: you
already did it once in a previous exercise.</p>
</div>
<p class="flip6">Click to Show/Hide Solution</p> <div class="panel6"><div class="highlight-python"><div class="highlight"><pre><span class="n">spectrum</span> <span class="o">=</span> <span class="n">img_bkg</span><span class="p">[</span><span class="mi">250</span><span class="p">:</span><span class="mi">260</span><span class="p">,</span> <span class="p">:]</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">spectrum</span><span class="p">)</span>
</pre></div>
</div>
</div><p><strong>To do</strong>: flux calibration and wavelength calibration!</p>
</div>
</div>
<div class="section" id="scipy">
<h1>SciPy<a class="headerlink" href="#scipy" title="Permalink to this headline">¶</a></h1>
<p>It is impossible to do justice to the full contents of the <a class="reference internal" href="#scipy">SciPy</a> package: is
entirely too large! What is left as homework for the reader is to
click through to the main <a class="reference external" href="http://docs.scipy.org/doc/scipy/reference/">SciPy Reference Manual</a> and skim the <a class="reference external" href="http://docs.scipy.org/doc/scipy/reference/tutorial/index.html">tutorial</a>. Keep
this repository of functionality in mind whenever you need some numerical
functionality that isn’t in NumPy: there is a good chance it is in SciPy:</p>
<ul class="simple">
<li>Basic functions in Numpy (and top-level scipy)</li>
<li>Special functions (scipy.special)</li>
<li>Integration (scipy.integrate)</li>
<li>Optimization (optimize)</li>
<li>Interpolation (scipy.interpolate)</li>
<li>Fourier Transforms (scipy.fftpack)</li>
<li>Signal Processing (signal)</li>
<li>Linear Algebra</li>
<li>Statistics</li>
<li>Multi-dimensional image processing (ndimage)</li>
<li>File IO (scipy.io)</li>
<li>Weave</li>
</ul>
</div>
</div>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
>index</a></li>
<li class="right" >
<a href="../contest/bounce.html" title="CONTEST: Make a fun bouncing balls demo"
>next</a> |</li>
<li class="right" >
<a href="ipython.html" title="IPython"
>previous</a> |</li>
<li><a href="../index.html">Python4Astronomers 1.1 documentation</a> »</li>
<li><a href="core.html" >Core packages for analysis: IPython, NumPy, and SciPy</a> »</li>
</ul>
</div>
<small>Copyright:
Smithsonian Astrophysical Observatory under terms of
<a rel="license"
href="http://creativecommons.org/licenses/by/3.0/">CC
Attribution 3.0</a></small> <a rel="license"
href="http://creativecommons.org/licenses/by/3.0/"><img alt="Creative Commons
License" style="border-width:0"
src="http://i.creativecommons.org/l/by/3.0/80x15.png" /></a>
</body>
</html>