diff --git a/src/NASABrightnessTemperature.jl b/src/NASABrightnessTemperature.jl index 7c1b6ac..18b8dd3 100644 --- a/src/NASABrightnessTemperature.jl +++ b/src/NASABrightnessTemperature.jl @@ -1,5 +1,40 @@ module NASABrightnessTemperature -# Write your package code here. +## Modules Used +using Logging +using NetRC +using Printf +using Statistics +import Base: download, show, read + +## Reexporting exported functions within these modules +using Reexport +@reexport using Dates +@reexport using GeoRegions +@reexport using NCDatasets + +## Exporting the following functions: +export + download, read, setup, extract, smoothing + + +modulelog() = "$(now()) - NASABrightnessTemperature.jl" + +function __init__() + setup() +end + +## Including Relevant Files + +include("setup.jl") + +include("dataset.jl") +include("download.jl") +include("save.jl") +include("read.jl") +include("extract.jl") +include("smoothing.jl") +include("filesystem.jl") +include("backend.jl") end diff --git a/src/backend.jl b/src/backend.jl new file mode 100644 index 0000000..c0b0235 --- /dev/null +++ b/src/backend.jl @@ -0,0 +1,93 @@ +## DateString Aliasing +yrmo2dir(date::TimeType) = Dates.format(date,dateformat"yyyy/mm") +yrmo2str(date::TimeType) = Dates.format(date,dateformat"yyyymm") +yr2str(date::TimeType) = Dates.format(date,dateformat"yyyy") +ymd2str(date::TimeType) = Dates.format(date,dateformat"yyyymmdd") + +function btdlonlat() + lon = convert(Array,-179.95:0.1:179.95) + lat = convert(Array,-89.95:0.1:89.95) + return lon,lat +end + +function checkdates( + dtbeg :: TimeType, + dtend :: TimeType +) + + if dtbeg < Date(2000,6,1) + error("$(modulelog()) - You have specified a date that is before the earliest available date of GPM IMERG data, 2000-06-01.") + end + + if dtend > (Dates.now() - Day(3)) + error("$(modulelog()) - You have specified a date that is later than the latest available date of GPM IMERG Near-Realtime data, $(now() - Day(3)).") + end + +end + +function ncoffsetscale(data::AbstractArray{<:Real}) + + init = data[findfirst(!isnan,data)] + dmax = init + dmin = init + for ii in eachindex(data) + dataii = data[ii] + if !isnan(dataii) + if dataii > dmax; dmax = dataii end + if dataii < dmin; dmin = dataii end + end + end + + scale = (dmax-dmin) / 65533; + offset = (dmax+dmin-scale) / 2; + + return scale,offset + +end + +function real2int16!( + outarray :: AbstractArray{Int16}, + inarray :: AbstractArray{<:Real}, + scale :: Real, + offset :: Real +) + + for ii in eachindex(inarray) + + idata = (inarray[ii] - offset) / scale + if isnan(idata) + outarray[ii] = -32767 + else; outarray[ii] = round(Int16,idata) + end + + end + + return + +end + +function nanmean( + data :: AbstractArray, + dNaN :: AbstractArray +) + nNaN = length(dNaN) + for iNaN in 1 : nNaN + dNaN[iNaN] = !isnan(data[iNaN]) + end + dataii = @view data[dNaN] + if !isempty(dataii); return mean(dataii); else; return NaN; end +end + +function nanmean( + data :: AbstractArray, + dNaN :: AbstractArray, + wgts :: AbstractArray, +) + nNaN = length(dNaN) + for iNaN in 1 : nNaN + dNaN[iNaN] = !isnan(data[iNaN]) + end + dataii = view(data,dNaN) .* view(wgts,dNaN) + wgtsii = view(wgts,dNaN) + if !isempty(dataii); return sum(dataii) / sum(wgtsii); else; return NaN; end +end \ No newline at end of file diff --git a/src/dataset.jl b/src/dataset.jl new file mode 100644 index 0000000..2a9fa8f --- /dev/null +++ b/src/dataset.jl @@ -0,0 +1,59 @@ +struct TbDataset{ST<:AbstractString, DT<:TimeType} + ID :: ST + name :: ST + doi :: ST + start :: DT + stop :: DT + path :: ST + hroot :: ST + fpref :: ST + fsuff :: ST +end + +function TbDataset( + ST = String, + DT = Date; + start :: TimeType = Date(2000,6), + stop :: TimeType = Dates.now() - Day(2), + path :: AbstractString = homedir(), +) + + @info "$(modulelog()) - Setting up data structure containing information on Early IMERG Daily data to be downloaded" + + fol = joinpath(path,"imergearlydy"); if !isdir(fol); mkpath(fol) end + fol = joinpath(path,"imergmask"); if !isdir(fol); mkpath(fol) end + + start = Date(year(start),month(start),1) + stop = Date(year(stop),month(stop),daysinmonth(stop)) + checkdates(start,stop) + + return TbDataset{ST,DT}( + "tb", "Brightness Temperature", "10.5067/P4HZB9N27EKU", + start, stop, + joinpath(path,"tb"), + "https://disc2.gesdisc.eosdis.nasa.gov/opendap/MERGED_IR/GPM_MERGIR.1", + "merg", "4km-pixel.nc4", + ) + +end + +function show( + io :: IO, + btd :: TbDataset{ST,DT} +) where {ST<:AbstractString, DT<:TimeType} + + print( + io, + "The NASA Brightness Temperature Dataset {$ST,$DT} has the following properties:\n", + " Dataset ID (ID) : ", btd.ID, '\n', + " Logging Name (name) : ", btd.name, '\n', + " DOI URL (doi) : ", btd.doi, '\n', + " Data Directory (path) : ", btd.path, '\n', + " Date Begin (start) : ", btd.start, '\n', + " Date End (stop) : ", btd.stop , '\n', + " Timestep : Half Hourly\n", + " Data Resolution : 4 km\n", + " Data Server (hroot) : ", btd.hroot, '\n', + ) + +end \ No newline at end of file diff --git a/src/download.jl b/src/download.jl new file mode 100644 index 0000000..080c391 --- /dev/null +++ b/src/download.jl @@ -0,0 +1,118 @@ +""" + download( + btd :: NASAPrecipitationDataset, + geo :: GeoRegion = GeoRegion("GLB"); + overwrite :: Bool = false + ) -> nothing + +Downloads a NASA Precipitation dataset specified by `btd` for a GeoRegion specified by `geo` + +Arguments +========= +- `btd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range +- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat + +Keyword Arguments +================= +- `overwrite` : If `true`, overwrite any existing files +""" +function download( + btd :: TbDataset{ST,DT}, + geo :: GeoRegion = GeoRegion("GLB"); + overwrite :: Bool = false +) where {ST<:AbstractString, DT<:TimeType} + + @info "$(modulelog()) - Downloading $(btd.name) data for the $(geo.name) GeoRegion from $(btd.start) to $(btd.stop)" + + lon,lat = btdlonlat(); nlon = length(lon) + ginfo = RegionGrid(geo,lon,lat) + + @info "$(modulelog()) - Preallocating temporary arrays for extraction of $(btd.name) data for the $(geo.name) GeoRegion from the original gridded dataset" + glon = ginfo.lon; nglon = length(glon); iglon = ginfo.ilon + glat = ginfo.lat; nglat = length(glat); iglat = ginfo.ilat + tmp0 = zeros(Float32,nglat,nglon) + var = zeros(Float32,nglon,nglat,48) + + if typeof(geo) <: PolyRegion + msk = ginfo.mask + else; msk = ones(nglon,nglat) + end + + if iglon[1] > iglon[end] + shift = true + iglon1 = iglon[1] : nlon; niglon1 = length(iglon1) + iglon2 = 1 : iglon[end]; niglon2 = length(iglon2) + tmp1 = @view tmp0[:,1:niglon1] + tmp2 = @view tmp0[:,niglon1.+(1:niglon2)] + @info "Temporary array sizes: $(size(tmp1)), $(size(tmp2))" + else + shift = false + iglon = iglon[1] : iglon[end] + end + + if iglat[1] > iglat[end] + iglat = iglat[1] : -1 : iglat[end] + else + iglat = iglat[1] : iglat[end] + end + + if btd.v6; varID = "precipitationCal"; else; varID = "precipitation" end + + for dt in btd.start : Day(1) : btd.stop + + fnc = btdfnc(btd,geo,dt) + if overwrite || !isfile(fnc) + + @info "$(modulelog()) - Downloading $(btd.name) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(dt) ..." + + ymdfnc = Dates.format(dt,dateformat"yyyymmdd") + btddir = joinpath(btd.hroot,"$(year(dt))",@sprintf("%03d",dayofyear(dt))) + + for it = 1 : 48 + + @debug "$(modulelog()) - Loading data into temporary array for timestep $(dyfnc[it])" + + btdfnc = "$(btd.fpref).$ymdfnc-$(dyfnc[it]).$(btd.fsuff)" + + tryretrieve = 0 + ds = 0 + while !(typeof(ds) <: NCDataset) && (tryretrieve < 20) + if tryretrieve > 0 + @info "$(modulelog()) - Attempting to request data from NASA OPeNDAP servers on Attempt $(tryretrieve+1) of 20" + end + ds = NCDataset(joinpath(btddir,btdfnc)) + tryretrieve += 1 + end + + if !shift + NCDatasets.load!(ds[varID].var,tmp0,iglat,iglon,1) + else + NCDatasets.load!(ds[varID].var,tmp1,iglat,iglon1,1) + NCDatasets.load!(ds[varID].var,tmp2,iglat,iglon2,1) + end + close(ds) + + @debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion" + for ilat = 1 : nglat, ilon = 1 : nglon + varii = tmp0[ilat,ilon] + mskii = msk[ilon,ilat] + if (varii != -9999.9f0) && !isnan(mskii) + var[ilon,ilat,it] = varii / 3600 + else; var[ilon,ilat,it] = NaN32 + end + end + end + + save(var,dt,btd,geo,ginfo) + + else + + @info "$(modulelog()) - $(btd.name) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(dt) exists in $(fnc), and we are not overwriting, skipping to next timestep ..." + + end + + flush(stderr) + + end + +end \ No newline at end of file diff --git a/src/filesystem.jl b/src/filesystem.jl new file mode 100644 index 0000000..a3b8a8a --- /dev/null +++ b/src/filesystem.jl @@ -0,0 +1,60 @@ +""" + npdfnc( + npd :: NASAPrecipitationDataset, + geo :: GeoRegion, + dt :: TimeType + ) -> String + +Returns of the path of the file for the NASA Precipitation dataset specified by `npd` for a GeoRegion specified by `geo` at a date specified by `dt`. + +Arguments +========= +- `npd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range +- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat +- `dt` : A specified date. The NCDataset retrieved may will contain data for the date, although it may also contain data for other dates depending on the `NASAPrecipitationDataset` specified by `npd` +""" +function btdfnc( + btd :: TbDataset{ST,DT}, + geo :: GeoRegion, + dt :: TimeType +) where {ST<:AbstractString, DT<:TimeType} + + fol = joinpath(btd.path,geo.ID,yrmo2dir(dt)) + fnc = btd.ID * "-" * geo.ID * "-" * ymd2str(dt) * ".nc" + return joinpath(fol,fnc) + +end + +#### + +function btdanc( + btd :: TbDataset{ST,DT}, + geo :: GeoRegion, + dt :: TimeType +) where {ST<:AbstractString, DT<:TimeType} + + fol = joinpath(btd.path,geo.ID) + fnc = btd.ID * "-" * geo.ID * "-" * yr2str(dt) * ".nc" + return joinpath(fol,fnc) + +end + +#### + +function btdsmth( + btd :: TbDataset{ST,DT}, + geo :: GeoRegion, + dt :: TimeType, + smoothlon :: Real, + smoothlat :: Real, + smoothtime :: Int +) where {ST<:AbstractString, DT<:TimeType} + + fol = joinpath(btd.path,geo.ID,yrmo2dir(dt)) + fnc = btd.ID * "-" * geo.ID * "-" * "smooth" * "_" * + @sprintf("%.2f",smoothlon) * "x" * @sprintf("%.2f",smoothlat) * + "_" * @sprintf("%02d",smoothtime) * "steps" * + "-" * ymd2str(dt) * ".nc" + return joinpath(fol,fnc) + +end \ No newline at end of file diff --git a/src/read.jl b/src/read.jl new file mode 100644 index 0000000..77cb539 --- /dev/null +++ b/src/read.jl @@ -0,0 +1,69 @@ +""" + read( + btd :: NASAPrecipitationDataset, + geo :: GeoRegion, + dt :: TimeType; + lonlat :: Bool = false + ) -> NCDataset (if lonlat = false) + -> lon, lat, NCDataset (if lonlat = true) + +Reads a NASA Precipitation dataset specified by `btd` for a GeoRegion specified by `geo` at a date specified by `dt`. + +Arguments +========= +- `btd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range +- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat +- `dt` : A specified date. The NCDataset retrieved may will contain data for the date, although it may also contain data for other dates depending on the `NASAPrecipitationDataset` specified by `btd` + +Keyword Arguments +================= +- `lonlat` : if `true`, then return the longitude and latitude vectors for the dataset. Otherwise only the NCDataset type will be returned. +""" +function read( + btd :: TbDataset, + geo :: GeoRegion, + dt :: TimeType; + smooth :: Bool = false, + smoothlon :: Real = 0, + smoothlat :: Real = 0, + smoothtime :: Real = 0, + quiet :: Bool = false +) + + bnc = btdfnc(btd,geo,dt) + + raw = true + if smooth + if iszero(smoothlon) && iszero(smoothlat) && iszero(smoothtime) + error("$(modulelog()) - Incomplete specification of smoothing parameters in either the longitude or latitude directions") + end + bnc = btdsmth(btd,geo,dt,smoothlon,smoothlat,smoothtime) + raw = false + end + + if quiet + disable_logging(Logging.Warn) + end + + if raw + if !isfile(bnc) + error("$(modulelog()) - The $(btd.name) Dataset for the $(geo.ID) GeoRegion at Date $dt does not exist at $(bnc). Check if files exist at $(btd.datapath) or download the files here") + end + @info "$(modulelog()) - Opening the $(btd.name) NCDataset in the $(geo.ID) GeoRegion for $dt" + end + if smooth + if !isfile(bnc) + error("$(modulelog()) - The spatially smoothed ($(@sprintf("%.2f",smoothlon))x$(@sprintf("%.2f",smoothlat))) $(btd.name) Dataset for $(geo.ID) GeoRegion at Date $dt does not exist at $(bnc). Check if files exist at $(btd.datapath) or download the files here") + end + @info "$(modulelog()) - Opening the spatialtemporally smoothed ($(@sprintf("%.2f",smoothlon))ºx$(@sprintf("%.2f",smoothlat))º, $(@sprintf("%02d",smoothtime)) timesteps) $(btd.name) NCDataset in the $(geo.ID) GeoRegion for $dt" + end + + if quiet + disable_logging(Logging.Debug) + end + + flush(stderr) + + return NCDataset(bnc) + +end \ No newline at end of file diff --git a/src/setup.jl b/src/setup.jl new file mode 100644 index 0000000..6d38963 --- /dev/null +++ b/src/setup.jl @@ -0,0 +1,68 @@ +""" + setup(; + login :: AbstractString = "", + password :: AbstractString = "", + overwrite :: Bool = false + ) + +Function that sets up your .dodsrc and .netrc files for you. If you do not provide the keyword arguments for `login` and `password`, the .netrc file will not be created. + +Keyword arguments +================= +- `login` : Your Earthdata username. You need to specify both `login` and `password` arguments. +- `password` : Your Earthdata password. You need to specify both `login` and `password` arguments. +- `overwrite` : If `true`, overwrite existing `.dodsrc` and `.netrc` files with the data provided. +""" +function setup(; + login :: AbstractString = "", + password :: AbstractString = "", + overwrite :: Bool = false +) + + fdodsrc = joinpath(homedir(),".dodsrc") + if !isfile(fdodsrc) + @info "$(modulelog()) - Setting up .dodsrc file for NASA OPeNDAP servers to point at cookie and .netrc directories ..." + open(fdodsrc,"w") do f + write(f,"HTTP.COOKIEJAR=$(joinpath(homedir(),".urs_cookies"))\n") + write(f,"HTTP.NETRC=$(joinpath(homedir(),".netrc"))") + end + else + if overwrite + @warn "$(modulelog()) - .dodsrc file exists at $fdodsrc, overwriting again" + open(fdodsrc,"w") do f + write(f,"HTTP.COOKIEJAR=$(joinpath(homedir(),".urs_cookies"))\n") + write(f,"HTTP.NETRC=$(joinpath(homedir(),".netrc"))") + end + else + @info "$(modulelog()) - .dodsrc file exists at $fdodsrc" + end + end + + if !netrc_check() + if login == "" && password == "" + @warn "$(modulelog()) - .netrc file does not exist at $(netrc_file()), you need to setup the .netrc file in order for NASAPrecipitation.jl to work" + end + end + + if netrc_check() && !netrc_checkmachine(netrc_read(logging=false),machine="urs.earthdata.nasa.gov") + if login == "" && password == "" + @warn "$(modulelog()) - No existing machine urs.earthdata.nasa.gov in .netrc file, please setup login and password for this machine for NASAPrecipitation.jl to work" + end + end + + if login != "" && password != "" + if overwrite || !netrc_check() + @info "$(modulelog()) - Setting up .netrc file containing login and password information for NASA OPeNDAP servers ..." + netrc = netrc_read() + if netrc_check() && netrc_checkmachine(netrc,machine="urs.earthdata.nasa.gov") + netrc_modify!(netrc,machine="urs.earthdata.nasa.gov",login=login,password=password) + else + netrc_add!(netrc,machine="urs.earthdata.nasa.gov",login=login,password=password) + end + netrc_write(netrc) + else + @info "$(modulelog()) - Existing .netrc file detected, since overwrite options is not selected, leaving file be ..." + end + end + +end \ No newline at end of file