|
| 1 | +if occursin(r"RvSpectMLEcoSystem$", pwd()) |
| 2 | + cd("NeidSolarScripts") |
| 3 | + using Pkg |
| 4 | + Pkg.activate(".") |
| 5 | + end |
| 6 | + |
| 7 | +verbose = true |
| 8 | + make_plots = true |
| 9 | + if verbose && !isdefined(Main,:RvSpectML) println("# Loading RvSpecML") end |
| 10 | + using RvSpectMLBase |
| 11 | + using EchelleInstruments, EchelleInstruments.NEID |
| 12 | + using EchelleCCFs |
| 13 | + #= |
| 14 | + using RvSpectML |
| 15 | + using NeidSolarScripts |
| 16 | + using NeidSolarScripts.SolarRotation |
| 17 | + if verbose println("# Loading other packages") end |
| 18 | + =# |
| 19 | + using CSV, DataFrames, Query, StatsBase, Statistics, Dates |
| 20 | + using JLD2, FileIO |
| 21 | + using NaNMath |
| 22 | + |
| 23 | +target_subdir = "good_days/DRPv0.7" # USER: Replace with directory of your choice |
| 24 | + fits_target_str = "Sun" |
| 25 | + output_dir = "output/continuum" |
| 26 | + #outputs = Dict{String,Any}() |
| 27 | + paths_to_search_for_param = [pwd(),"examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/scratch/jpn23/"] |
| 28 | + # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below |
| 29 | + pipeline_plan = PipelinePlan() |
| 30 | + dont_make_plot!(pipeline_plan, :movie) |
| 31 | + |
| 32 | +reset_all_needs!(pipeline_plan) |
| 33 | +#if need_to(pipeline_plan,:read_spectra) |
| 34 | +if verbose println("# Finding what data files are avaliable.") end |
| 35 | +if isfile("manifest.csv") |
| 36 | + df_files = CSV.read("manifest.csv", DataFrame) |
| 37 | + @assert size(df_files,1) >= 1 |
| 38 | + @assert hasproperty(df_files,:Filename) |
| 39 | + @assert hasproperty(df_files,:target) |
| 40 | + @assert hasproperty(df_files,:bjd) |
| 41 | + @assert hasproperty(df_files,:ssbz) |
| 42 | + @assert hasproperty(df_files,:exptime) |
| 43 | +else |
| 44 | + eval(read_data_paths(paths_to_search=paths_to_search_for_param)) |
| 45 | + @assert isdefined(Main,:neid_data_path) |
| 46 | + df_files = make_manifest(neid_data_path, target_subdir, NEID ) |
| 47 | + CSV.write("manifest.csv", df_files) |
| 48 | +end |
| 49 | + |
| 50 | +if verbose println("# Reading in customized parameters from param.jl.") end |
| 51 | + if !@isdefined(idx_day_to_use) |
| 52 | + idx_day_to_use = 1 |
| 53 | + end |
| 54 | + eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param)) |
| 55 | + if match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] == match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1] |
| 56 | + date_str = match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] |
| 57 | + else |
| 58 | + date_str = string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]) * "-" * string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1]) |
| 59 | + end |
| 60 | + #= |
| 61 | + outputs["df_files_use"] = df_files_use |
| 62 | +
|
| 63 | + outputs_filename = joinpath(output_dir,"solar_" * date_str * "_new.jld2") |
| 64 | + if isfile(outputs_filename) && false |
| 65 | + times_already_processed = load(outputs_filename, "times") |
| 66 | + files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1) |
| 67 | + if files_in_day_to_process == length(times_already_processed) |
| 68 | + println("# Already processed all ", length(times_already_processed), " files for ", date_str) |
| 69 | + exit() |
| 70 | + end |
| 71 | + end |
| 72 | + =# |
| 73 | + |
| 74 | +using Distributed |
| 75 | +addprocs(4) |
| 76 | + |
| 77 | +@everywhere using RCall |
| 78 | +@everywhere afs_src = joinpath(pwd(),"src","AFS.R") |
| 79 | +@everywhere R"source($afs_src)" |
| 80 | +@everywhere using Pkg |
| 81 | +@everywhere Pkg.activate(".") |
| 82 | +@everywhere using RvSpectMLBase |
| 83 | +@everywhere using EchelleInstruments |
| 84 | + |
| 85 | +@everywhere function calc_continuum_model(spectrum::AbstractSpectra2D; order_idx::Integer ) |
| 86 | + possible_pix = get_pixel_range(get_inst(spectrum),order_idx) |
| 87 | + bad_pix = bad_col_ranges(get_inst(spectrum),order_idx) |
| 88 | + pix_rng = EchelleInstruments.calc_complement_index_ranges(possible_pix,bad_pix) |
| 89 | + pix = mapreduce(p->collect(p),vcat,pix_rng) |
| 90 | + afs_inputs = zeros(Float64,length(pix),2) |
| 91 | + afs_inputs[:,1] .= spectrum.λ[pix,order_idx] |
| 92 | + afs_inputs[:,2] .= spectrum.flux[pix,order_idx] |
| 93 | + @assert !any(isnan.(afs_inputs)) |
| 94 | + #= |
| 95 | + wv = mapreduce(p->spec.λ[p,order_idx],vcat,pix_rng) |
| 96 | + @assert !any(isnan.(wv)) |
| 97 | + inten = mapreduce(p->convert(Vector{Float64},spec.flux[p,order_idx]),vcat,pix_rng) |
| 98 | + @assert !any(isnan.(inten)) |
| 99 | + afs_inputs = hcat(wv,inten) |
| 100 | + =# |
| 101 | + #df = DataFrame("wv"=>wv,"intes"=>inten) |
| 102 | + afs_output_R = R"AFS($afs_inputs,0.95,0.25)" |
| 103 | + afs_output = rcopy(afs_output_R) |
| 104 | + continuum = zeros(eltype(spectrum.flux),size(spectrum.flux[:,order_idx])) |
| 105 | + continuum = fill(NaN,size(spectrum.flux[:,order_idx])) |
| 106 | + continuum[pix] .= afs_output |
| 107 | + return continuum |
| 108 | +end |
| 109 | + |
| 110 | + |
| 111 | + |
| 112 | +@everywhere using EchelleInstruments.NEID |
| 113 | +@everywhere function calc_continuum_model(spectrum::AbstractSpectra2D ) |
| 114 | + vec_of_orders = pmap(ord->calc_continuum_model(spectrum,order_idx=ord), min_order(get_inst(spectrum)):max_order(get_inst(spectrum)) ) |
| 115 | + output = fill(NaN, size(spectrum.flux)) |
| 116 | + for (i,ord) in enumerate(min_order(get_inst(spectrum)):max_order(get_inst(spectrum))) |
| 117 | + output[:,ord] .= vec_of_orders[i] |
| 118 | + end |
| 119 | + return output |
| 120 | +end |
| 121 | + |
| 122 | +num_days_to_process = size(df_files_solar_by_day,1) |
| 123 | +for idx_day_to_use in 1:num_days_to_process |
| 124 | + df_files_use = df_files_solar_by_day[idx_day_to_use,:data] |> @orderby(_.bjd) |> DataFrame |
| 125 | + |
| 126 | + println("# *** Working on day ", idx_day_to_use, " with ", size(df_files_use,1), "files. ***" ) |
| 127 | + for row in eachrow(df_files_use) |
| 128 | + spec = NEID.read_data(row) |
| 129 | + m = match(r"neidL1_(\d+)[T_](\d+)\.fits$",row.Filename) |
| 130 | + output_filename = "neidL1_" * m.captures[1] * "T" * m.captures[2] * ".jld2" |
| 131 | + output_filename = joinpath(output_dir,output_filename) |
| 132 | + println("# Working on ", output_filename) |
| 133 | + continuum = calc_continuum_model(spec) |
| 134 | + @save output_filename continuum |
| 135 | + end |
| 136 | + |
| 137 | +end |
0 commit comments