Merge branch 'CyanideCN:master' into mergeCR

.github/workflows/ci.yml

@@ -31,7 +31,7 @@ jobs:
         run: |
           python -m pip install --upgrade setuptools
           python -m pip install --upgrade wheel
-          python -m pip install --upgrade numpy==2.0.0rc1
+          python -m pip install --upgrade numpy>=2.0
           python setup.py bdist_wheel
       - name: Upload artifact
         uses: actions/upload-artifact@v2

.github/workflows/test.yml

@@ -0,0 +1,31 @@
+name: Test
+
+on:
+  push:
+    branches: [ master ]
+  pull_request:
+    branches: [ master ]
+
+jobs:
+  test:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+    - uses: conda-incubator/setup-miniconda@v2
+      with:
+        auto-update-conda: true
+        python-version: "3.11"
+        environment-file: environment.yaml
+        activate-environment: build
+    - name: Install
+      run: |
+        python -m pip install --upgrade setuptools
+        python -m pip install --upgrade numpy>=2.0
+        python -m pip install .
+    - name: Run tests
+      run: |
+        python -m pip install pytest
+        pytest

cinrad/io/_dtype.py

@@ -11,7 +11,8 @@
     'S_SPECIAL_dtype', 'CC2_header', 'CC2_obs', 'CC2_data', 'CC2_other', 'PA_radial',
     'L3_raster','L3_hail','L3_meso','L3_feature','L3_tvs','L3_sti_header','L3_sti_motion',
     'L3_sti_position','L3_sti_attribute','L3_sti_component','L3_sti_adaptation',
-    'L3_vwp_header','L3_vwp','L3_swp','L3_uam','mocm_dtype','mocm_si_dtype','mocm_si_block'
+    'L3_vwp_header','L3_vwp','L3_swp','L3_uam','mocm_dtype','mocm_si_dtype','mocm_si_block',
+    'L3_wer_header'
 ]
 # fmt: on
 from cinrad.io._radar_struct.CC import (
@@ -44,6 +45,7 @@
     l3_vwp_table_dtype as L3_vwp,
     l3_swp as L3_swp,
     l3_uam_dtype as L3_uam,
+    l3_wer_header_dtype as L3_wer_header,
 )
 from cinrad.io._radar_struct.CD import (
     radarsite_dtype as CD_site,

cinrad/io/_radar_struct/standard_data.py

@@ -399,3 +399,12 @@
 ]
 
 l3_uam_dtype = np.dtype(l3_uam)
+
+l3_wer_header = [
+    ("elevation", "f4"),
+    ("scan_time", "i4"),
+    ("center_height", "i4"),
+    ("res", "20c"),
+]
+
+l3_wer_header_dtype = np.dtype(l3_wer_header)

cinrad/io/level3.py

@@ -196,6 +196,7 @@ class SWAN(object):
     Args:
         file (str, IO): Path points to the file or a file object.
     """
+
     dtype_conv = {0: "B", 1: "b", 2: "u2", 3: "i2", 4: "u2"}
     size_conv = {0: 1, 1: 1, 2: 2, 3: 2, 4: 2}
 
@@ -484,12 +485,14 @@ def __init__(self, file):
         if self.name == "None":
             self.name = self.code
         del self.geo
-        if self.ptype in [1, 13, 14, 25, 26, 27, 28, 51, 52]:  # PPI radial format
+        if self.ptype in [1, 13, 14, 24, 25, 26, 27, 28, 51, 52]:  # PPI radial format
             self._parse_radial_fmt()
-        elif self.ptype in [6, 8, 9, 10, 18, 23]:
+        elif self.ptype in [4, 6, 8, 9, 10, 18, 23]:
             self._parse_raster_fmt()
         elif self.ptype == 3:
             self._parse_cappi_fmt()
+        elif self.ptype == 20:
+            self._parse_wer_fmt()
         elif self.ptype == 32:
             self._parse_vwp_fmt()
         elif self.ptype == 36:
@@ -504,6 +507,10 @@ def __init__(self, file):
             self._parse_tvs_fmt()
         elif self.ptype == 44:
             self._parse_uam_fmt()
+        else:
+            raise RadarDecodeError(
+                "Unsupported product type {}:{}".format(self.ptype, self.pname)
+            )
         self.f.close()
 
     def _parse_header(self):
@@ -523,7 +530,9 @@ def _parse_header(self):
         ph = np.frombuffer(self.f.read(128), SDD_pheader)
         self.ptype = ph["product_type"][0]
         self.pname = self.ptype_corr[self.ptype]
-        self.scantime = datetime.datetime.utcfromtimestamp(ph["scan_start_time"][0])
+        self.scantime = datetime.datetime.fromtimestamp(
+            ph["scan_start_time"][0], datetime.timezone.utc
+        )
         if self.ptype == 1:
             self.pname = self.dtype_corr[ph["dtype_1"][0]]
         self.params = ProductParamsParser.parse(self.ptype, self.f.read(64))
@@ -698,18 +707,40 @@ def _parse_cappi_fmt(self):
         self._dataset = ds
 
     def _parse_vwp_fmt(self):
-        vwp_header = np.frombuffer(self.f.read(32), L3_vwp_header)
-        vwp_count = vwp_header["number_of_vols"][0]
-        vwps = np.frombuffer(self.f.read(32 * vwp_count), L3_vwp)
-        height = DataArray(vwps["height"])
-        wind_direction = DataArray(vwps["wind_direction"])
-        wind_speed = DataArray(vwps["wind_speed"])
+        self.vwp_header = np.frombuffer(self.f.read(32), L3_vwp_header)
+        timestamp = list()
+        height = list()
+        wd = list()
+        ws = list()
+        rms = list()
+        while True:
+            buf = self.f.read(32)
+            if not buf:
+                break
+            vwp = np.frombuffer(buf, L3_vwp)
+            timestamp.append(vwp["start_time"][0])
+            height.append(vwp["height"][0])
+            wd.append(vwp["wind_direction"][0])
+            ws.append(vwp["wind_speed"][0])
+            rms.append(vwp["rms_std"][0])
+        height = list(set(height))
+        timestamp = list(set(timestamp))
+        height.sort()
+        timestamp.sort()
+        shape = (len(timestamp), len(height))
+        wd = np.round(np.array(wd).astype(float).reshape(shape), 0)
+        ws = np.round(np.array(ws).astype(float).reshape(shape), 2)
+        rms = np.round(np.array(rms).astype(float).reshape(shape), 2)
+        wd_da = DataArray(
+            wd,
+            coords=[
+                timestamp,
+                height,
+            ],
+            dims=["times", "height"],
+        )
         ds = Dataset(
-            {
-                "wind_direction": wind_direction,
-                "wind_speed": wind_speed,
-                "height": height,
-            },
+            {"wind_direction": wd_da},
             attrs={
                 "scan_time": self.scantime.strftime("%Y-%m-%d %H:%M:%S"),
                 "site_code": self.code,
@@ -719,6 +750,8 @@ def _parse_vwp_fmt(self):
                 "task": self.task_name,
             },
         )
+        ds["wind_speed"] = (["times", "height"], ws)
+        ds["rms"] = (["times", "height"], rms)
         self._dataset = ds
 
     def _parse_swp_fmt(self):
@@ -1002,6 +1035,23 @@ def _parse_uam_fmt(self):
         ds["latitude"] = DataArray(lat[:, 0])
         self._dataset = ds
 
+    def _parse_wer_fmt(self):
+        wer = Dataset()
+        while True:
+            buf = self.f.read(32)
+            if not buf:
+                break
+            wer_header = np.frombuffer(buf, L3_wer_header)
+            elev = wer_header["elevation"][0]
+            self._parse_raster_fmt()
+            if len(wer) == 0:
+                wer = self._dataset.copy()
+                wer = wer.rename({self.pname: "{}_{:.1f}".format(self.pname, elev)})
+                wer.attrs["center_height"] = wer_header["center_height"][0]
+            else:
+                wer["{}_{:.1f}".format(self.pname, elev)] = self._dataset[self.pname]
+        self._dataset = wer
+
     def get_data(self):
         return self._dataset
 

cinrad/visualize/utils.py

@@ -17,7 +17,9 @@
 from cartopy.io import shapereader
 import cartopy.crs as ccrs
 from cartopy.mpl.geoaxes import GeoAxes
+from cartopy.feature import Feature
 import shapefile
+import shapely.geometry as sgeom
 from vanadis.colormap import Colormap
 from cinrad_data import get_font_path, get_shp_list, get_shp_file
 
@@ -138,9 +140,9 @@ def _get_uniform_cmap(cmap: Any) -> Any:
                      "0.9", "0.7", "0.5", "0.3", "0.1"],
              "OHP":["", "203.2", "152.4", "101.6", "76.2", "63.5", "50.8", "44.45", "38.1", "31.75",
                     "25.4", "19.05", "12.7", "6.35", "2.54", "0"],
-             "cHCL":["Drizzle", "Rain", "Ice Crystals", "Dry Snow", "Wet Snow", "Vertical Ice", 
+             "cHCL":["Drizzle", "Rain", "Ice Crystals", "Dry Snow", "Wet Snow", "Vertical Ice",
                      "Low-Dens Graupel", "High-Dens Graupel", "Hail", "Big Drops", ""],
-             "HCL":["Rain", "Heavy Rain", "Hail", "Big Drops", "Clear-Air Echo", "Ground Clutter", 
+             "HCL":["Rain", "Heavy Rain", "Hail", "Big Drops", "Clear-Air Echo", "Ground Clutter",
                     "Dry snow", "Wet snow", "Ice Crystals", "Graupel", "Unknown", ""]}
 # fmt: on
 
@@ -272,6 +274,19 @@ def get_shp() -> list:
     return shps
 
 
+class _ShapelyFeature(Feature):
+    r"""Copied from cartopy.feature.ShapelyFeature"""
+
+    def __init__(self, geometries, crs, **kwargs):
+        super().__init__(crs, **kwargs)
+        if isinstance(geometries, sgeom.base.BaseGeometry):
+            geometries = [geometries]
+        self._geoms = tuple(geometries)
+
+    def geometries(self):
+        return iter(self._geoms)
+
+
 def add_shp(
     ax: Any,
     proj: ccrs.Projection,
@@ -287,29 +302,35 @@ def add_shp(
         line_colors = ["grey", "lightgrey", "white"]
     elif style == "white":
         line_colors = ["lightgrey", "grey", "black"]
-    ax.add_geometries(
-        shps[0],
-        shp_crs,
-        edgecolor=line_colors[0],
-        facecolor="None",
-        zorder=3,
-        linewidth=0.5,
+    ax.add_feature(
+        _ShapelyFeature(
+            shps[0],
+            shp_crs,
+            edgecolor=line_colors[0],
+            facecolor="None",
+            zorder=3,
+            linewidth=0.5,
+        )
     )
-    ax.add_geometries(
-        shps[1],
-        shp_crs,
-        edgecolor=line_colors[1],
-        facecolor="None",
-        zorder=3,
-        linewidth=0.7,
+    ax.add_feature(
+        _ShapelyFeature(
+            shps[1],
+            shp_crs,
+            edgecolor=line_colors[1],
+            facecolor="None",
+            zorder=3,
+            linewidth=0.7,
+        )
     )
-    ax.add_geometries(
-        shps[2],
-        shp_crs,
-        edgecolor=line_colors[2],
-        facecolor="None",
-        zorder=3,
-        linewidth=1,
+    ax.add_feature(
+        _ShapelyFeature(
+            shps[2],
+            shp_crs,
+            edgecolor=line_colors[2],
+            facecolor="None",
+            zorder=3,
+            linewidth=1,
+        )
     )
     if coastline:
         ax.coastlines(resolution="10m", color=line_colors[2], zorder=3, linewidth=1)

setup.py

@@ -51,11 +51,11 @@
         "vanadis",
         "cinrad_data>=0.1"
     ],
-    package_data={
-        data_pth: glob.glob(join(data_pth, "*.*")),
-        join(data_pth, "colormap"): glob.glob(join(data_pth, "colormap", "*.cmap")),
-        join("cinrad", "correct"): [join("cinrad", "correct", "unwrap_2d_ljmu.c")],
-    },
+    package_dir={"cinrad": "cinrad"},
+    package_data={"cinrad": [
+        "data/*.*",
+        "data/*/*.*"
+    ]},
     scripts=[],
     ext_modules=ext_modules,
     include_dirs=include_dirs,