Crude support for CfRadial

README

@@ -6,7 +6,7 @@
 
 --------------------------------------------------
 
-Versioon 2023.3.30
+Versioon 2023.5.17
 
 Tarmo Tanilsoo, 2023
 
@@ -58,6 +58,7 @@ Vaata ka demovideot, kus näidatakse töö käiku: http://youtu.be/RZOYDvJv434
 * Python Imaging Library, võib olla ka selle Pillow fork(arendusmasinal jooksebki juba see, kuna see tuleb uuemate Ubuntu distributsioonidega kaasa). Muuhulgas peab olema ka ImageTk moodul, vähemalt Ubuntulistel on see eraldi pakis.
 * H5Py
 * numpy
+* netCDF4
 
 Soovitatavalt vähemalt 1 GB RAM'i.
 
@@ -74,6 +75,7 @@ Windowsi kasutajatel peaks olema ilmselt kõige lihtsam paigaldada käsurealt Py
 * IRIS RAW - samad produktid, mis ODIM H5 puhul. (fail peab lõppema .raw laiendiga!)
 * DORADE - väga toores tugi.
 * JMA GRIB - polaarkoordinaadistikus andmed. Faili nimi peab järgima JMA struktuuri, kuna seda kasutatakse produkti tüübi tuvastamiseks.
+* CfRadial - veel toores, sihitud PyART'i poolt genereeritud failidele
 
 Olen avatud täiendavatele formaatidele toe lisamisele, kuid sellisel juhul vajan formaadi dokumentatsiooni, näidisfaile ning näidet oodatavast väljundist.
 

README_ENGLISH

@@ -6,7 +6,7 @@
 
 --------------------------------------------------
 
-Version 2023.3.30
+Version 2023.5.17
 
 Tarmo Tanilsoo, 2023
 _____
@@ -59,8 +59,9 @@ as it is distributed with newer Ubuntu distributions). Also needs to have ImageT
 package at least on Ubuntu)
 * H5PY
 * Numpy
+* NetCDF4
 
-If you're using Windows and Tkinter is installed, then the easiest way is probably to fetch packages h5py and Pillow from Python Package Index repository. Other dependencies should be coming automatically with them.
+If you're using Windows and Tkinter is installed, then the easiest way is probably to fetch packages h5py and Pillow from Python Package Index repository. Other dependencies should be mostly coming automatically with them.
 
 I recommend to have at least 1 GB of RAM.
 
@@ -77,6 +78,7 @@ I recommend to have at least 1 GB of RAM.
 * IRIS RAW - file name must end with .raw extension
 * DORADE - very crude support at the moment.
 * JMA GRIB - polar coordinate data. File name must follow JMA's convention for correct product detection.
+* CfRadial - still quite crude, aimed at files generated by PyART
 
 I am open to adding additional formats, but for that I'd need the format documentation and example files to experiment with, as well as an example of expected output.
 

src/decoderadar.py

@@ -44,7 +44,7 @@
 import sys
 import platform
 from copy import deepcopy
-
+import netCDF4
 from array import array
 import numpy as np
 from numpy import fromstring
@@ -57,6 +57,93 @@
 NP_INT32 = np.int32
 NP_FLOAT = np.float32
 
+class CfRadial():
+    def processTime(self,timestamp): #Convert CfRadial timestamp to Python datetime
+        return datetime.datetime.strptime(timestamp.decode("utf-8").strip().replace("\x00",""), "%Y-%m-%dT%H:%M:%SZ")
+
+    def __init__(self,path):
+        self.type = "CfRadial"
+        self.headers = {}
+        self.azimuths = []
+        self.data = []
+        self.times = []
+        self.elevations = []
+        self.elevationNumbers = []
+        self.quantities = []
+        self.isModified = False
+        self._raw = netCDF4.Dataset(path)
+
+        self._raw.set_auto_maskandscale(False)
+
+        self.nominalElevations=np.round_(self._raw["fixed_angle"][:],2).tolist()
+        if "time_reference" in self._raw.variables:
+            referenceTime = self.processTime(self._raw["time_reference"][:].tobytes())
+        else:
+            referenceTime = self.processTime(self._raw["time_coverage_start"][:].tobytes())
+
+        #Wavelength calculation
+        prtratio = self._raw["prt_ratio"][:][0]
+        prt = self._raw["prt"][:][0]
+        highprf = round(1/prt,2)
+        lowprf = round(highprf/prtratio,2)
+        rmax = self._raw["unambiguous_range"][:][0]
+        vmax = self._raw["nyquist_velocity"][:][0]
+        for numerator in range(1,10):
+            if round(prtratio*numerator,2) % 1 == 0: break
+        self.wavelength = ((vmax/numerator)*4)/highprf
+
+        self.headers={"timestamp":referenceTime,
+                      "latitude":float(self._raw.variables["latitude"][0]),
+                      "longitude":float(self._raw.variables["longitude"][0]),
+                      "height": float(self._raw.variables["altitude"][0])}
+
+        sweepStarts = self._raw.variables["sweep_start_ray_index"][:].tolist()
+        sweepEnds = self._raw.variables["sweep_end_ray_index"][:].tolist()
+        nSweeps = len(sweepStarts)
+
+        #Conversion between radial quantity names and ODIM
+
+        quantityNames = {
+            "total_power": "TH",
+            "reflectivity": "DBZH",
+            "velocity": "VRADH",
+            "spectrum_width": "WRADH",
+            "differential_reflectivity": "ZDR",
+            "specific_differential_phase": "KDP",
+            "differential_phase": "PHIDP",
+            "normalized_coherent_power":"NCP",
+            "cross_correlation_ratio":"RHOHV",
+            "radar_echo_classification":"HCLASS",
+            "corrected_velocity":"VRADDH"}
+
+        for i in range(len(sweepStarts)):
+            times = []
+            for t in self._raw["time"][sweepStarts[i]:sweepEnds[i]+1][:]:
+                times.append(referenceTime + datetime.timedelta(seconds = t))
+            self.times.append(times)
+            self.elevations.append(self._raw["elevation"][sweepStarts[i]:sweepEnds[i]+1][:])
+            self.azimuths.append(self._raw["azimuth"][sweepStarts[i]:sweepEnds[i]+1][:])
+            dataEntry = {}
+            quantitiesInSweep = []
+            for j in quantityNames:
+                quantityCode = quantityNames[j]
+                if j in self._raw.variables:
+                    quantitiesInSweep.append(quantityCode)
+                    sweepData = np.array(self._raw[j][sweepStarts[i]:sweepEnds[i]+1])
+                    dataEntry[quantityCode] = {"data" : sweepData,
+                                               "dataType": np.uint16,
+                                               "rscale": self._raw["range"].meters_between_gates / 1000,
+                                               "rstart": (self._raw["range"].meters_to_center_of_first_gate - (self._raw["range"].meters_between_gates / 2)) / 1000, #Coordinate system used by TRV assumes the rstart indicates the start of bin, not centre of bin
+                                               "highprf":highprf,
+                                               "lowprf":lowprf,
+                                               "offset":None,
+                                               "gain":None,
+                                               "nodata":0,
+                                               "undetect":self._raw[j]._FillValue}
+            self.data.append(dataEntry)
+            self.quantities.append(quantitiesInSweep)
+            self.elevationNumbers.append(len(self.elevations))
+
 
 def JMARLE(data, NBITS, MAXV):
     LNGU = 2**(NBITS)-1-MAXV
@@ -1625,10 +1712,6 @@ def __init__(self,path):
         del(sisu) #Garbage collecting - this file can be bigly!
 
 class HDF5():
-    def warn(self,message):
-        if message not in self.issuedWarnings:
-            print(message,file=sys.stderr)
-            self.issuedWarnings.append(message)
     def __init__(self,path):
         self.type="HDF5"
         self.headers={}
@@ -1639,7 +1722,6 @@ def __init__(self,path):
         self.elevations=[]
         self.quantities=[]
         self.isModified = False
-        self.issuedWarnings=[]
 
         andmed=HDF5Fail(path,"r") #Load the data file
 
@@ -2679,7 +2761,10 @@ def leiasuund(rad,rad2,y,currentDisplay,zoom=1,centre=[1000,1000],samm=1):
 
 def scaleValue(value,gain,offset,nodata,undetect,rangefolding):
     if value != nodata and value != undetect and value != rangefolding:
-        return round(value*gain+offset,6)
+        if gain and offset:
+            return round(value*gain+offset,6)
+        else: #In case we have already floats (like in CfRadial)
+            return value
     else:
         if value == rangefolding:
             return "RF"

src/main.pyw

@@ -541,11 +541,16 @@ class BatchExportWindow(Tkinter.Toplevel): #Window for batch export
                 path=self.datadir+"/"+i
                 currentfilepath=path
                 stream=file_read(path)
-                if path[-3:]== ".h5" or stream[1:4]==b"HDF":
+                if path[-3:]== ".h5":
                     productChoice.config(state=Tkinter.NORMAL)
                     elevationChoice.config(state=Tkinter.NORMAL)
                     hcanames=fraasid["iris_hca"]
                     currentlyOpenData=HDF5(path)
+                elif path[-3:]== ".nc":
+                    productChoice.config(state=Tkinter.NORMAL)
+                    elevationChoice.config(state=Tkinter.NORMAL)
+                    hcanames=fraasid["iris_hca"]
+                    currentlyOpenData=CfRadial(path)
                 elif stream[0:4] == b"AR2V" or stream[0:8] == b"ARCHIVE2":
                     productChoice.config(state=Tkinter.NORMAL)
                     elevationChoice.config(state=Tkinter.NORMAL)
@@ -1831,11 +1836,16 @@ def load(path=None,defaultElevation=0):
             productChoice.config(state=Tkinter.NORMAL)
             elevationChoice.config(state=Tkinter.NORMAL)
             currentlyOpenData=JMA(path)
-        elif path[-3:]== ".h5" or stream[1:4]==b"HDF":
+        elif path[-3:]== ".h5":
             productChoice.config(state=Tkinter.NORMAL)
             elevationChoice.config(state=Tkinter.NORMAL)
             hcanames=fraasid["iris_hca"]
             currentlyOpenData=HDF5(path)
+        elif path[-3:]== ".nc":
+            productChoice.config(state=Tkinter.NORMAL)
+            elevationChoice.config(state=Tkinter.NORMAL)
+            hcanames=fraasid["iris_hca"]
+            currentlyOpenData=CfRadial(path)
         elif path[-4:].lower() == ".raw":
             productChoice.config(state=Tkinter.NORMAL)
             elevationChoice.config(state=Tkinter.NORMAL)

src/translations.py

@@ -145,7 +145,7 @@ def fixArabic(nas="مرحبا"):
 phrases={"estonian":
          {
              "LANG_ID":"EE",
-             "name":"TRV 2023.3.30",
+             "name":"TRV 2023.5.17",
              "loading_states":"Laen andmeid... osariigid",
              "coastlines":"Rannajooned",
              "countries":"Maismaapiirid",
@@ -327,7 +327,7 @@ def fixArabic(nas="مرحبا"):
          "english":
          {
              "LANG_ID":"EN",
-             "name":"TRV 2023.3.30",
+             "name":"TRV 2023.5.17",
              "loading_states":"Loading data... states",
              "coastlines":"Coastlines",
              "countries":"Countries",
@@ -509,7 +509,7 @@ def fixArabic(nas="مرحبا"):
          "arabic":  #NOTE: I am not a native speaker. Corrections are more than welcome as pull requests, especially from radar meteorologists(even from PME!)
          {
              "LANG_ID":"AR",
-             "name":u"TRV 2023.3.30",
+             "name":u"TRV 2023.5.17",
              "loading_states":u"فتح البيانات غيوغرافي... ولايات",
              "coastlines":u"خطوط الساحل",
              "countries":u"دول",
@@ -692,7 +692,7 @@ def fixArabic(nas="مرحبا"):
          "japanese":
          {
              "LANG_ID":"JP",
-             "name":"TRV 2023.3.30",
+             "name":"TRV 2023.5.17",
              "loading_states":u"データをロード中... 州",
              "coastlines":u"海岸線",
              "countries":u"郡の境界",