添加了部分StandardPUP产品解析

cinrad/io/_dtype.py

@@ -9,7 +9,7 @@
     'SDD_cut', 'SDD_rad_header', 'SDD_mom_header', 'SAB_dtype', 'CAB_dtype',
     'SWAN_dtype', 'CD_dtype', 'CD_DATA', 'SDD_pheader', 'L3_radial', 'L3_rblock',
     'S_SPECIAL_dtype', 'CC2_header', 'CC2_obs', 'CC2_data', 'CC2_other', 'PA_radial',
-    'L3_raster'
+    'L3_raster','HAIL_dtype'
 ]
 # fmt: on
 from cinrad.io._radar_struct.CC import (
@@ -135,3 +135,14 @@
         ("rec", _CD_record, 998),
     ]
 )
+
+PUP_HAIL_TABLE=[
+    ("hail_id", "i4"),
+    ("hail_azimuth", "f4"),
+    ("hail_range", "i4"),
+    ("hail_possibility", "i4"),
+    ("hail_severe_possibility", "i4"),
+    ("hail_size", "f4"),
+    ("rcm", "i4"),
+]
+HAIL_dtype= np.dtype(PUP_HAIL_TABLE)

cinrad/io/level3.py

@@ -430,7 +430,9 @@ def _empty(buf):
     def parse(cls, product_type, param_bytes):
         buf = BytesIO(param_bytes)
         map_func = {1: cls._ppi, 2: cls._rhi, 51: cls._ppi, 18: cls._empty}
-        params = map_func[product_type](buf)
+        params = {"elevation": 0}
+        if product_type in map_func.keys():
+            params = map_func[product_type](buf)
         buf.close()
         return params
 
@@ -442,6 +444,7 @@ class StandardPUP(RadarBase):
                   9:'RHO', 10:'PHI', 11:'KDP', 12:'CP', 14:'HCL', 15:'CF', 16:'SNRH',
                   17:'SNRV', 32:'Zc', 33:'Vc', 34:'Wc', 35:'ZDRc', 71:'RR', 72:'HGT',
                   73:'VIL', 74:'SHR', 75:'RAIN', 76:'RMS', 77:'CTR'}
+    ptype_corr = {1:'REF', 6:'ET', 18:'CR', 25:'OHP', 26:'THP', 27:'STP', 28:'USP', 38:'HI'}
     # fmt: on
     def __init__(self, file):
         self.f = prepare_file(file)
@@ -450,10 +453,12 @@ def __init__(self, file):
         self.stationlat = self.geo["lat"][0]
         self.stationlon = self.geo["lon"][0]
         self.radarheight = self.geo["height"][0]
-        if self.ptype == 1:  # PPI radial format
+        if self.ptype in [1, 25, 26, 27, 28]:  # PPI radial format
             self._parse_radial_fmt()
-        elif self.ptype == 18:
+        elif self.ptype in [6, 18]:
             self._parse_raster_fmt()
+        elif self.ptype == 38:
+            self._parse_hail_fmt()
         if self.name == "None":
             self.name = self.code
         del self.geo
@@ -475,6 +480,7 @@ def _parse_header(self):
         self.scan_config = np.frombuffer(self.f.read(256 * cut_num), SDD_cut)
         ph = np.frombuffer(self.f.read(128), SDD_pheader)
         self.ptype = ph["product_type"][0]
+        self.ptype_name = self.ptype_corr[self.ptype]
         self.scantime = datetime.datetime.utcfromtimestamp(ph["scan_start_time"][0])
         self.dtype = self.dtype_corr[ph["dtype_1"][0]]
         if self.dtype in ["Zc", "REF"] and self.ptype == 18:
@@ -489,9 +495,10 @@ def _parse_radial_fmt(self):
         reso = radial_header["reso"][0] / 1000
         start_range = radial_header["start_range"][0] / 1000
         end_range = radial_header["max_range"][0] / 1000
+        nradial = radial_header["nradial"][0]
         data = list()
         azi = list()
-        while True:
+        for i in range(nradial):
             buf = self.f.read(32)
             if not buf:
                 break
@@ -507,6 +514,7 @@ def _parse_radial_fmt(self):
         data_rf = np.ma.masked_not_equal(raw, 1)
         raw = np.ma.masked_less(raw, 5)
         data = (raw - offset) / scale
+        data = np.ma.masked_less(data, 1)
         az = np.linspace(0, 360, raw.shape[0])
         az += azi[0]
         az[az > 360] -= 360
@@ -522,7 +530,7 @@ def _parse_radial_fmt(self):
         )
         da = DataArray(data, coords=[azi, dist], dims=["azimuth", "distance"])
         ds = Dataset(
-            {self.dtype: da},
+            {self.ptype_name: da},
             attrs={
                 "elevation": self.params["elevation"],
                 "range": end_range,
@@ -570,7 +578,7 @@ def _parse_raster_fmt(self):
             dims=["latitude", "longitude"],
         )
         ds = Dataset(
-            {self.dtype: da},
+            {self.ptype_name: da},
             attrs={
                 "elevation": 0,
                 "range": max_range,
@@ -583,6 +591,48 @@ def _parse_raster_fmt(self):
             },
         )
         self._dataset = ds
-
+
+        def _parse_hail_fmt(self):
+            hail_count = np.frombuffer(self.f.read(4), "i4")[0]
+            hail_table = np.frombuffer(self.f.read(hail_count * 28), HAIL_dtype)
+            ht0msl = np.frombuffer(self.f.read(4), "f4")[0]
+            ht20msl = np.frombuffer(self.f.read(4), "f4")[0]
+            hail_azimuth = hail_table["hail_azimuth"]
+            hail_range = hail_table["hail_range"]
+            hail_size = DataArray(hail_table["hail_size"])
+            hail_possibility = DataArray(hail_table["hail_possibility"])
+            hail_severe_possibility = DataArray(hail_table["hail_severe_possibility"])
+            lon, lat = list(), list()
+            for dis, az in zip(hail_range, hail_azimuth):
+                x, y = get_coordinate(
+                    dis / 1000,
+                    az * deg2rad,
+                    self.params["elevation"],
+                    self.stationlon,
+                    self.stationlat,
+                )
+                lon.append(x)
+                lat.append(y)
+            ds = Dataset(
+                {
+                    "hail_possibility": hail_possibility,
+                    "hail_size": hail_size,
+                    "hail_severe_possibility": hail_severe_possibility,
+                },
+                attrs={
+                    "scan_time": self.scantime.strftime("%Y-%m-%d %H:%M:%S"),
+                    "site_code": self.code,
+                    "site_name": self.name,
+                    "site_longitude": self.stationlon,
+                    "site_latitude": self.stationlat,
+                    "task": self.task_name,
+                    "height_0deg": ht0msl,
+                    "height_-20deg": ht20msl,
+                },
+            )
+            ds["longitude"] = DataArray(lon)
+            ds["latitude"] = DataArray(lat)
+            self._dataset = ds
+
     def get_data(self):
         return self._dataset