-
Notifications
You must be signed in to change notification settings - Fork 8
/
io.jl
178 lines (161 loc) · 5.23 KB
/
io.jl
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
export load, save
"""
load(fitsfile::String, IntensityMap)
This loads in a fits file that is more robust to the various imaging algorithms
in the EHT, i.e. is works with clean, smili, eht-imaging.
The function returns an tuple with an intensitymap and a second named tuple with ancillary
information about the image, like the source name, location, mjd, and radio frequency.
"""
function fileio_load(file::File{format"FITS"}, T::Type{IntensityMap})
return _load_fits(file.filename, T)
end
function _load_fits(fname, ::Type{IntensityMap})
img, head = FITS(fname, "r") do f
if ndims(f[1])==2
image = copy(read(f[1])')[:,end:-1:1]
else
@warn "Currently only stokes I is loaded"
image = Matrix{Float64}(read(f[1])[:,:,1,1]')
end
header = read_header(f[1])
nx = Int(header["NAXIS1"])
ny = Int(header["NAXIS2"])
psizex = abs(float(header["CDELT1"]))*π/180
psizey = abs(float(header["CDELT2"]))*π/180
ra = float(header["OBSRA"])
dec = float(header["OBSDEC"])
#Get frequency
freq = 0.0
if haskey(header, "FREQ")
freq = parse(Float64, string(header["FREQ"]))
elseif "CRVAL3" in keys(header)
freq = float(header["CRVAL3"])
end
mjd = 0.0
if haskey(header, "MJD")
mjd = parse(Float64, string(header["MJD"]))
end
source = "NA"
if haskey(header,"OBJECT")
source = string(header["OBJECT"])
end
bmaj = 1.0 #Nominal values
bmin = 1.0
if haskey(header, "BUNIT")
if header["BUNIT"] == "JY/BEAM"
@info "Converting Jy/Beam => Jy/pixel"
bmaj = header["BMAJ"]*π/180
bmin = header["BMIN"]*π/180
beamarea = (2.0*π*bmaj*bmin)/(8*log(2))
image .= image.*(psizex*psizey/beamarea)
end
end
imap = IntensityMap(image, psizex*nx, psizey*ny)
info = (source=source, RA=ra, DEC=dec, mjd=mjd, ν=freq)
return imap, info
end
end
"""
save(file::String, img::IntensityMap, obs)
Saves an image to a fits file. You can optionally pass an EHTObservation so that ancillary information
will be added.
"""
function fileio_save(fname::File{format"FITS"}, img::IntensityMap, obs = nothing)
head = make_header(obs)
_save_fits(fname.filename, img, head)
end
function make_header(obs)
if isnothing(obs)
return (source="NA", RA=0.0, DEC=0.0, mjd=0, freq=0.0)
else
return (source=String(obs.source), RA=obs.ra, DEC=obs.dec, mjd=obs.mjd, freq=obs.frequency)
end
end
function _save_fits(fname::String, image::IntensityMap, head)
headerkeys = ["SIMPLE",
"BITPIX",
"NAXIS",
"NAXIS1",
"NAXIS2",
"OBJECT",
"CTYPE1",
"CTYPE2",
"CDELT1",
"CDELT2",
"OBSRA",
"OBSDEC",
"FREQ",
"CRPIX1",
"CRPIX2",
"MJD",
"TELESCOP",
"BUNIT",
"STOKES"]
psizex, psizey = pixelsizes(image)
values = [true,
-64,
2,
size(image, 1),
size(image, 2),
head.source,
"RA---SIN",
"DEC---SIN",
rad2deg(psizex),
rad2deg(psizey),
head.RA,
head.DEC,
head.freq,
size(image,1)/2+0.5,
size(image,2)/2+0.5,
head.mjd,
"VLBI",
"JY/PIXEL",
"STOKES"]
comments = ["conforms to FITS standard",
"array data type",
"number of array dimensions",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
"",
""]
FITS(fname, "w") do hdu
hdeheader = FITSHeader(headerkeys, values, comments)
img = copy(image[:,end:-1:1]')
write(hdu, img, header=hdeheader)
end
end
"""
$(SIGNATURES)
Load a ThemisPy style ascii EHT observation file.
"""
function load_tpy(file)
data = readdlm(file, skipstart=1)
bs1 = Symbol.(getindex.(data[:,5], Ref(1:2)))
bs2 = Symbol.(getindex.(data[:,5], Ref(3:4)))
baselines = tuple.(bs1, bs2)
edata = StructArray{EHTVisibilityDatum{Float64}}(
visr=float.(data[:,8]),
visi=float.(data[:,10]),
error=float.(data[:,9]),
u=float.(data[:,6])*1e6,
v=float.(data[:,7])*1e6,
time=float.(data[:,4]),
frequency=fill(227e9, size(data,1)),
bandwidth=fill(4e6, size(data,1)),
baseline=baselines
)
mjd = Int(modified_julian(UTCEpoch(Int(data[1,2]), 1,1,0)).Δt+float(data[1,3]))
return EHTObservation(data=edata, mjd=mjd, ra=180.0, dec=0.0, source=Symbol(data[1,1]))
end