Adding some basic stuff, based off NASAPrecipitation.jl

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

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

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

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

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