Managing FITS files is a big deal if you're an observer. Is the interface in Julia any good for this?
Task:
minimum functionThen we'll show how you can easily compare the speeds of operations in Julia, and use this to check whether it's faster to do this by using Astropy or using FITSIO.
The fits file used was downloaded from NASA's sample FITS page (the HST WFPC II file).
In [1]:
# set up plotting
using PyPlot
cmap = get_cmap("gist_stern")
INFO: Loading help data...
Out[1]:
In [2]:
# load in astropy
using PyCall
@pyimport astropy.io.fits as pyf
In [3]:
# load in the data
hdulist = pyf.open("data/fits-image.fits")
Out[3]:
1-element Array{Any,1}:
PyObject <astropy.io.fits.hdu.image.PrimaryHDU object at 0x7f0e7ac99450>
In [4]:
# first and only HDU
hdu = hdulist[1]
header = hdu["header"]
Out[4]:
PyObject SIMPLE = T / file does conform to FITS standard
BITPIX = -32 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 100 / length of data axis 1
NAXIS2 = 100 / length of data axis 2
EXTEND = T / FITS dataset may contain extensions
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
ORIGIN = 'NOAO-IRAF FITS Image Kernel December 2001' / FITS file originator
DATE = '2003-07-24T18:28:58' / Date FITS file was generated
IRAF-TLM= '10:28:58 (24/07/2003)' / Time of last modification
OBJECT = 'U5780205B[1/4]' / Name of the object observed
CRVAL1 = 182.62570646887
CRVAL2 = 39.41233768771
CRPIX1 = -2.135000000000000E+02
CRPIX2 = -2.040000000000000E+02
CD1_1 = 2.756050E-5
CD1_2 = -2.082210E-6
CD2_1 = -2.080210E-6
CD2_2 = -2.758710E-5
MIR_REVR= T
ORIENTAT= -174.798
FILLCNT = 0
ERRCNT = 0
FPKTTIME= 51229.799079
LPKTTIME= 51229.799246
CTYPE1 = 'RA---TAN'
CTYPE2 = 'DEC--TAN'
DETECTOR= 4
DEZERO = 311.7072
BIASEVEN= 311.7648
BIASODD = 311.6496
GOODMIN = -3.289304
GOODMAX = 3420.24
DATAMEAN= 0.7084697
GPIXELS = 628289
SOFTERRS= 0
CALIBDEF= 1816
STATICD = 0
ATODSAT = 57
DATALOST= 0
BADPIXEL= 0
OVERLAP = 0
PHOTMODE= 'WFPC2,4,A2D7,FR533P15,,CAL '
PHOTFLAM= 8.566962E-19
PHOTZPT = -21.1
PHOTPLAM= 5782.972
PHOTBW = 2070.218
MEDIAN = -0.06821804
MEDSHADO= -0.03253492
HISTWIDE= 0.9845759
SKEWNESS= -4.585337
MEANC10 = 0.1359666
MEANC25 = 0.04015671
MEANC50 = 0.1129934
MEANC100= 1.897491
MEANC200= 0.9687142
MEANC300= 1.05377
BACKGRND= 0.3467301
FILETYPE= 'SCI ' / type of data found in data file
TELESCOP= 'HST' / telescope used to acquire data
INSTRUME= 'WFPC2 ' / identifier for instrument used to acquire data
EQUINOX = 2000.0 / equinox of celestial coord. system
/ WFPC-II DATA DESCRIPTOR KEYWORDS
ROOTNAME= 'u5780205r ' / rootname of the observation set
PROCTIME= 5.284443195602E+04 / Pipeline processing time (MJD)
OPUS_VER= 'OPUS 14.3b ' / OPUS software system version number
CAL_VER = ' ' / CALWP2 code version
/ SCIENCE INSTRUMENT CONFIGURATION
MODE = 'FULL' / instr. mode: FULL (full res.), AREA (area int.)
SERIALS = 'OFF' / serial clocks: ON, OFF
/ IMAGE TYPE CHARACTERISTICS
IMAGETYP= 'EXT ' / DARK/BIAS/IFLAT/UFLAT/VFLAT/KSPOT/EXT/ECAL
CDBSFILE= 'NO ' / GENERIC/BIAS/DARK/PREF/FLAT/MASK/ATOD/NO
PKTFMT = 96 / packet format code
/ FILTER CONFIGURATION
FILTNAM1= 'FR533P15' / first filter name
FILTNAM2= ' ' / second filter name
FILTER1 = 69 / first filter number (0-48)
FILTER2 = 0 / second filter number (0-48)
FILTROT = 15.0 / partial filter rotation angle (degrees)
LRFWAVE = 4877.000000 / linear ramp filter wavelength
/ INSTRUMENT STATUS USED IN DATA PROCESSING
UCH1CJTM= -88.2569 / TEC cold junction #1 temperature (Celsius)
UCH2CJTM= -88.6697 / TEC cold junction #2 temperature (Celsius)
UCH3CJTM= -88.3028 / TEC cold junction #3 temperature (Celsius)
UCH4CJTM= -88.7671 / TEC cold junction #4 temperature (Celsius)
UBAY3TMP= 13.2302 / bay 3 A1 temperature (deg C)
KSPOTS = 'OFF' / Status of Kelsall spot lamps: ON, OFF
SHUTTER = 'A' / Shutter in place at beginning of the exposure
ATODGAIN= 7.0 / Analog to Digital Gain (Electrons/DN)
/ RSDP CONTROL KEYWORDS
MASKCORR= 'COMPLETE' / Do mask correction: PERFORM, OMIT, COMPLETE
ATODCORR= 'COMPLETE' / Do A-to-D correction: PERFORM, OMIT, COMPLETE
BLEVCORR= 'COMPLETE' / Do bias level correction
BIASCORR= 'COMPLETE' / Do bias correction: PERFORM, OMIT, COMPLETE
DARKCORR= 'COMPLETE' / Do dark correction: PERFORM, OMIT, COMPLETE
FLATCORR= 'SKIPPED ' / Do flat field correction
SHADCORR= 'OMIT ' / Do shaded shutter correction
DOSATMAP= 'OMIT ' / Output saturated pixel map
DOPHOTOM= 'COMPLETE' / Fill photometry keywords
DOHISTOS= 'OMIT ' / Make histograms: PERFORM, OMIT, COMPLETE
OUTDTYPE= 'REAL ' / Output image datatype: REAL, LONG, SHORT
/ CALIBRATION REFERENCE FILES
MASKFILE= 'uref$f8213081u.r0h ' / name of the input DQF of known bad pixels
ATODFILE= 'uref$dbu1405iu.r1h' / name of the A-to-D conversion file
BLEVFILE= 'ucal$u5780205r.x0h ' / Engineering file with extended register da
BLEVDFIL= 'ucal$u5780205r.q1h ' / Engineering file DQF
BIASFILE= 'uref$j9a1612mu.r2h' / name of the bias frame reference file
BIASDFIL= 'uref$j9a1612mu.b2h' / name of the bias frame reference DQF
DARKFILE= 'uref$j2g1549cu.r3h' / name of the dark reference file
DARKDFIL= 'uref$j2g1549cu.b3h' / name of the dark reference DQF
FLATFILE= 'uref$f4i1559cu.r4h' / name of the flat field reference file
FLATDFIL= 'uref$f4i1559cu.b4h' / name of the flat field reference DQF
SHADFILE= 'uref$e371355eu.r5h' / name of the reference file for shutter sha
PHOTTAB = ' ' / name of the photometry calibration table
GRAPHTAB= 'mtab$n3d14531m_tmg.fits' / the HST graph table
COMPTAB = 'mtab$n431303km_tmc.fits' / the HST components table
/ DEFAULT KEYWORDS SET BY STSCI
SATURATE= 4095 / Data value at which saturation occurs
USCALE = 1.0 / Scale factor for output image
UZERO = 0.0 / Zero point for output image
/ READOUT DURATION INFORMATION
READTIME= 464 / Length of time for CCD readout in clock ticks
/ PLANETARY SCIENCE KEYWORDS
PA_V3 = 49.936909 / position angle of V3-axis of HST (deg)
RA_SUN = 3.337194516616E+02 / right ascension of the sun (deg)
DEC_SUN = -1.086675160382E+01 / declination of the sun (deg)
EQNX_SUN= 2000.0 / equinox of the sun
MTFLAG = F / moving target flag; T if it is a moving target
EQRADTRG= 0.000000 / equatorial radius of target (km)
FLATNTRG= 0.000000 / flattening of target
NPDECTRG= 0.000000 / north pole declination of target (deg)
NPRATRG = 0.000000 / north pole right ascension of target (deg)
ROTRTTRG= 0.000000 / rotation rate of target
LONGPMER= 0.000000 / longitude of prime meridian (deg)
EPLONGPM= 0.000000 / epoch of longitude of prime meridian (sec)
SURFLATD= 0.000000 / surface feature latitude (deg)
SURFLONG= 0.000000 / surface feature longitude (deg)
SURFALTD= 0.000000 / surface feature altitude (km)
/ PODPS FILL VALUES
PODPSFF = 0 / 0=(no podps fill); 1=(podps fill present)
STDCFFF = 0 / 0=(no st dcf fill); 1=(st dcf fill present)
STDCFFP = '0x5569' / st dcf fill pattern (hex)
RSDPFILL= -100 / bad data fill value for calibrated images
/ EXPOSURE TIME AND RELATED INFORMATION
UEXPODUR= 300 / commanded duration of exposure (sec)
NSHUTA17= 1 / Number of AP17 shutter B closes
DARKTIME= 3.000000000000E+02 / Dark time (seconds)
UEXPOTIM= 16880 / Major frame pulse time preceding exposure start
PSTRTIME= '1999.051:19:08:37 ' / predicted obs. start time (yyyy.ddd:hh:mm:ss)
PSTPTIME= '1999.051:19:16:37 ' / predicted obs. stop time (yyyy.ddd:hh:mm:ss)
/ EXPOSURE INFORMATION
SUNANGLE= 141.618347 / angle between sun and V1 axis
MOONANGL= 126.698997 / angle between moon and V1 axis
SUN_ALT = -31.523479 / altitude of the sun above Earth's limb
FGSLOCK = 'FINE ' / commanded FGS lock (FINE,COARSE,GYROS,UNKNOWN)
DATE-OBS= '1999-02-20' / UT date of start of observation (yyyy-mm-dd)
TIME-OBS= '19:03:13' / UT time of start of observation (hh:mm:ss)
EXPSTART= 5.122979390428E+04 / exposure start time (Modified Julian Date)
EXPEND = 5.122979737650E+04 / exposure end time (Modified Julian Date)
EXPTIME = 300. / exposure duration (seconds)--calculated
EXPFLAG = 'NORMAL ' / Exposure interruption indicator
/ TARGET & PROPOSAL ID
TARGNAME= 'NGC4151 ' / proposer's target name
RA_TARG = 1.826355000000E+02 / right ascension of the target (deg) (J2000)
DEC_TARG= 3.940576666667E+01 / declination of the target (deg) (J2000)
ECL_LONG= 164.096619 / ecliptic longitude of the target (deg) (J2000)
ECL_LAT = 36.623709 / ecliptic latitude of the target (deg) (J2000)
GAL_LONG= 155.079532 / galactic longitude of the target (deg) (J2000)
GAL_LAT = 75.062679 / galactic latitude of the target (deg) (J2000)
PROPOSID= 8019 / PEP proposal identifier
PEP_EXPO= '02-030 ' / PEP exposure identifier including sequence
LINENUM = '02.030 ' / PEP proposal line number
SEQLINE = ' ' / PEP line number of defined sequence
SEQNAME = ' ' / PEP define/use sequence name
HISTORY MASKFILE=uref$f8213081u.r0h MASKCORR=COMPLETED
HISTORY PEDIGREE=INFLIGHT 01/01/1994 - 15/05/1995
HISTORY DESCRIP=STATIC MASK - INCLUDES CHARGE TRANSFER TRAPS
HISTORY BIASFILE=uref$j9a1612mu.r2h BIASCORR=COMPLETED
HISTORY PEDIGREE=INFLIGHT 29/08/98 - 21/08/99
HISTORY DESCRIP=not significantly different from j6e16008u.
HISTORY DARKFILE=uref$j2g1549cu.r3h DARKCORR=COMPLETED
HISTORY PEDIGREE=INFLIGHT 16/02/1999 - 16/02/1999
HISTORY DESCRIP=Pipeline dark: 120 frame superdark with hotpixels from
HISTORY 16/02/99
HISTORY FLATFILE=uref$f4i1559cu.r4h FLATCORR=SKIPPED
HISTORY PEDIGREE=DUMMY 18/04/1995
HISTORY DESCRIP=All pixels set to value of 1. Not flat-fielded.
HISTORY PC1: bias jump level ~0.100 DN.
HISTORY The following throughput tables were used:
HISTORY crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY crwfpc2comp$wfpc2_dqepc1_005_syn.fits,
HISTORY crwfpc2comp$wfpc2_a2d7pc1_004_syn.fits,
HISTORY crwfpc2comp$wfpc2_flatpc1_003_syn.fits
HISTORY The following throughput tables were used:
HISTORY crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY crwfpc2comp$wfpc2_dqewfc2_005_syn.fits,
HISTORY crwfpc2comp$wfpc2_a2d7wf2_004_syn.fits,
HISTORY crwfpc2comp$wfpc2_flatwf2_003_syn.fits
HISTORY The following throughput tables were used:
HISTORY crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY crwfpc2comp$wfpc2_dqewfc3_005_syn.fits,
HISTORY crwfpc2comp$wfpc2_a2d7wf3_004_syn.fits,
HISTORY crwfpc2comp$wfpc2_flatwf3_003_syn.fits
HISTORY The following throughput tables were used:
HISTORY crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY crwfpc2comp$wfpc2_dqewfc4_005_syn.fits,
HISTORY crwfpc2comp$wfpc2_a2d7wf4_004_syn.fits,
HISTORY crwfpc2comp$wfpc2_flatwf4_003_syn.fits
HISTORY Ran task WARMPIX, version 1.1 (Sep 28, 1999) at Thu 10:24:40
HISTORY 24-Jul-2003, rej_thresh=0.10000, fix_thresh=0.00300,
HISTORY var_thresh=0.00300,
HISTORY fix_dqval=1024, rej_val=INDEF
SKYSUB1 = 0. / Sky value for chip 1
PYS_VERS= 1.8 / PyStack.py: Version
IRAF_I = 'IRAF V2.12.1 July 2002 release:2.12'
STAC_I = '$RCSfile: imstack_imshift_imcomb.cl,v $ $Revision: 1.1 $'
PYS_ITER= 5 / PyStack.py: Number of iterations
PYS_HEXP= 2 / PyStack.py: Value of exponent
SKYSUB2 = 0.01 / Sky value for chip 2
SKYSUB3 = -0.08 / Sky value for chip 3
SKYSUB4 = -0.03 / Sky value for chip 4
NCOMBINE= 2 / Number of exposures for this association
HISTORY Ran task WMOSAIC, version 2.1 (Jun 1995), at Thu 10:28:58 24-Jul-2003
COMMENT
COMMENT
In [5]:
imdata = hdu[:data]
Out[5]:
100x100 Array{Float32,2}:
1.03158 1.08554 0.511125 … 7.40503 6.80991 5.86466
0.476544 0.369005 0.524518 7.08262 7.12883 5.27363
1.26589 0.799219 0.218222 7.61593 9.10855 6.91347
1.18256 0.878865 0.204786 8.80005 8.40067 7.88117
0.927722 0.752777 0.584217 8.84981 8.04905 7.60208
0.482523 0.127957 0.27255 … 8.67788 7.67705 8.29861
0.406336 0.394538 0.268648 8.63436 6.78155 8.92089
0.940592 0.52642 0.691145 7.9926 6.79646 8.01739
1.27879 0.915581 0.52578 7.55549 6.94876 7.38695
0.697374 0.52315 0.314522 7.16162 6.69554 6.81852
-0.134354 0.515386 0.738036 … 7.27472 6.81387 6.00003
-0.0580475 0.665356 0.614793 7.2126 5.88157 4.65218
0.720517 0.68849 0.855741 8.33559 5.1659 4.01831
⋮ ⋱
1.00261 0.63835 0.377774 1.09655 0.759311 0.635015
0.547952 0.700689 0.269449 0.575319 -0.0219589 0.783337
0.66705 0.420551 0.0238154 … 0.684515 0.597013 0.83292
0.573887 0.39462 0.49563 0.600988 1.09629 0.58552
0.296847 0.0363498 0.250984 0.868481 0.936388 0.250051
0.50536 0.322373 -0.10824 0.83425 0.0869685 0.116338
-0.197808 -0.171628 -0.0166526 0.827145 0.898985 0.418266
0.279576 -0.739037 0.00439202 … 0.253484 0.553173 0.279392
1.05891 0.347345 -0.164634 0.251521 0.54144 0.456686
0.0482905 0.171759 -0.0425407 0.368141 0.481697 0.0305117
-0.0420618 0.59643 -0.387436 0.441201 0.498471 0.0837238
-0.0112627 -0.0142497 0.0181858 0.584082 0.412454 0.4871
In [6]:
imshow(imdata, cmap=cmap, interpolation="bicubic")
Out[6]:
PyObject <matplotlib.image.AxesImage object at 0x7f0e7ac462d0>
In [7]:
# get the minimum pixel value
minimum(imdata)
Out[7]:
-1.8990651f0
In [8]:
# get the position of the brightest pixel(s)
@show brightest = maximum(imdata)
find(x -> x==brightest, imdata)
brightest = maximum(imdata) => 1975.2683f0
Out[8]:
1-element Array{Int64,1}:
8310
In [9]:
# the find command gives single indexes: turn it into coords using ind2sub
ind2sub(size(imdata), 8310)
Out[9]:
(10,84)
In [10]:
# find all the negative pixels
negatives = find(x -> x<0, imdata)
# copy the array and set those pixels to 0
newimdata = copy(imdata)
newimdata[negatives] = 0
# write back to file
pyf.writeto("data/new-image.fits", newimdata, header, clobber=true)
WARNING: Overwriting existing file 'data/new-image.fits'. [astropy.io.fits.file]
In [11]:
using FITSIO
hdulist = FITS("data/fits-image.fits")
Out[11]:
file: data/fits-image.fits
mode: r
extnum exttype extname
1 image_hdu
In [12]:
image = hdulist[1]
Out[12]:
file: data/fits-image.fits
extension: 1
type: IMAGE
image info:
bitpix: -32
size: (100,100)
In [13]:
imdata = read(image)
imshow(imdata, cmap=cmap, interpolation="bicubic")
Out[13]:
PyObject <matplotlib.image.AxesImage object at 0x7f0e7a603910>
In [14]:
function read_data_with_python()
hdulist = pyf.open("data/fits-image.fits")
image = hdulist[1]
data = image[:data]
end
function read_data_with_fitsio()
hdulist = FITS("data/fits-image.fits")
image = hdulist[1]
data = read(image)
end
# run them once to compile before timing
read_data_with_python()
read_data_with_fitsio();
In [15]:
@time read_data_with_python()
@time read_data_with_fitsio();
elapsed time: 0.009074418 seconds (530944 bytes allocated)
elapsed time: 0.000551839 seconds (42008 bytes allocated)
Content source: swt30/ioa-julia-tutorials
Similar notebooks: