some updates, np.nan i.s.o. np.NaN, added coamps to spw but that will…

… be removed again at some point

cht_cyclones/__init__.py

@@ -6,3 +6,5 @@
 """
 
 __version__ = "0.0.1"
+
+from .tropical_cyclone import TropicalCyclone

cht_cyclones/fit_holland_2010.py

@@ -301,7 +301,7 @@ def compute_mean_error(r, w, obs, wrad):
             if not np.isnan(obs["quadrants_radii"][irad, iquad]):
                 err[iquad, irad] = vrad - wrad[irad]
             else:
-                err[iquad, irad] = np.NAN
+                err[iquad, irad] = np.nan
 
     # Get error values
     mask = ~np.isnan(err)  # Create the mask

cht_cyclones/spiderweb.py

@@ -19,8 +19,14 @@ def __init__(self):
         self.ds = xr.Dataset()
         self.ds.attrs["description"] = "Tropical cyclone spiderweb by cht_cyclones"
 
-    def read(self, filename, fmt):
-        pass
+    def read(self, filename):
+        # Read spw file
+        # Get file extension
+        fmt = filename.split(".")[-1]        
+        if fmt == "nc":
+            self.ds = xr.open_dataset(filename)
+        elif fmt == "spw":
+            self.read_spiderweb_ascii(filename)
 
     def initialize_grid(self, track, nrad, ndir, spiderweb_radius):
         # Make Xarray dataset
@@ -132,6 +138,127 @@ def write(self, filename, format="netcdf", background_pressure=1013.0, tref=None
                                        merge_frac=merge_frac,
                                        include_rainfall=include_rainfall)
 
+    def read_spiderweb_ascii(self, filename):
+        # Read in ASCII
+
+        # Open file
+        fid = open(filename, "r")
+
+        # First loop through all lines and count number of times that a line starts with "TIME"
+        n_times = 0
+        while True:
+            line = fid.readline()
+            if not line:
+                # End of file
+                break
+            if line.split()[0] == "n_rows":
+                n_rows = int(line.split()[2])
+            if line.split()[0] == "n_cols":
+                n_cols = int(line.split()[2])
+            if line.split()[0] == "n_quantity":
+                n_quantity = int(line.split()[2])    
+            if line.split()[0] == "TIME":
+                n_times += 1
+        # Rewind file
+        fid.seek(0)
+
+        # Initialize arrays
+        t = np.zeros(n_times, dtype="datetime64[s]")
+        lon = np.zeros(n_times)
+        lat = np.zeros(n_times)
+
+        # Now find the number of rows and columns
+        n_rows = 0
+        n_cols = 0
+        for i in range(10):
+            line = fid.readline()
+            if line.split()[0] == "n_rows":
+                n_rows = int(line.split()[2])
+            if line.split()[0] == "n_cols":
+                n_cols = int(line.split()[2])
+        # Rewind file
+        fid.seek(0)
+
+
+        # Read header information
+        vsn = fid.readline().split()[2]
+        filetype = fid.readline().split()[2]
+        nodata_value = float(fid.readline().split()[2])
+        n_cols = int(fid.readline().split()[2])
+        n_rows = int(fid.readline().split()[2])
+        grid_unit = fid.readline().split()[2]
+        spiderweb_radius = float(fid.readline().split()[2])
+        spw_rad_unit = fid.readline().split()[2]
+        spw_merge_frac = float(fid.readline().split()[2])
+        n_quantity = int(fid.readline().split()[2])
+        quantity1 = fid.readline().split()[2]
+        quantity2 = fid.readline().split()[2]
+        quantity3 = fid.readline().split()[2]
+        unit1 = fid.readline().split()[2]
+        unit2 = fid.readline().split()[2]
+        unit3 = fid.readline().split()[2]
+
+        if n_quantity == 4:
+            quantity4 = fid.readline().split()[2]
+            unit4 = fid.readline().split()[2]
+
+        # Initialize arrays
+        wind_speed = np.zeros((n_times, n_rows, n_cols))
+        wind_from_direction = np.zeros((n_times, n_rows, n_cols))
+        pressure_drop = np.zeros((n_times, n_rows, n_cols))
+        if n_quantity == 4:
+            precipitation = np.zeros((n_times, n_rows, n_cols))
+
+        for it in range(n_times):
+            line = fid.readline()
+            parts = line.split("=")
+            # Get the second part
+            line = parts[1]
+            t0 = float(line.split()[0])
+            tunits = line.split()[1]
+            trefstrs = line.split()[3:-1]
+            trefstr = " ".join(trefstrs)
+            tref = np.datetime64(trefstr)
+            if tunits[0:2] == "mi":
+                t = tref + np.timedelta64(int(t0), "m")
+            elif tunits[0:2] == "ho":
+                t = tref + np.timedelta64(int(t0), "h")
+            elif tunits[0:2] == "se":
+                t = tref + np.timedelta64(int(t0), "s")
+            x_spw_eye = float(fid.readline().split()[2])
+            y_spw_eye = float(fid.readline().split()[2])
+            p_drop_spw_eye = float(fid.readline().split()[2])
+
+
+            # Read the data
+            for i in range(n_rows):
+                wind_speed[it, i, :] = np.array(fid.readline().split(), dtype=float)
+            for i in range(n_rows):
+                wind_from_direction[it, i, :] = np.array(fid.readline().split(), dtype=float)
+            for i in range(n_rows):
+                pressure_drop[it, i, :] = np.array(fid.readline().split(), dtype=float)
+            if n_quantity == 4:
+                for i in range(n_rows):
+                    precipitation[it, i, :] = np.array(fid.readline().split(), dtype=float)
+
+        # Close file
+        fid.close()
+
+        # # Create Xarray dataset
+        # self.ds = xr.Dataset(
+        #     {
+        #         "wind_speed": (["time", "range", "azimuth"], wind_speed),
+        #         "wind_from_direction": (["time", "range", "azimuth"], wind_from_direction),
+        #         "pressure": (["time", "range", "azimuth"], 101200.0 - pressure_drop),
+        #     },
+        #     coords={
+        #         "time": np.arange(0, n_times),
+        #         "range": np.arange(0, n_rows),
+        #         "azimuth": np.arange(0, n_cols),
+        #     },
+        # )
+        # xxx=1
+
     def write_spiderweb_ascii(self, filename, background_pressure=1013.0, tref=None, merge_frac=0.5, include_rainfall=False):
         # Write in ASCII
 

cht_cyclones/track.py

@@ -88,94 +88,94 @@ def convert_units_imperial_to_metric(self, config):
                 if self.gdf.vmax[it] > 0.0:
                     self.gdf.vmax[it] = self.gdf.vmax[it] * knots_to_ms * config["wind_conversion_factor"]
                 else:
-                    self.gdf.vmax[it] = np.NaN    
+                    self.gdf.vmax[it] = np.nan    
                 if self.gdf.rmw[it] > 0.0:
                     self.gdf.rmw[it] = self.gdf.rmw[it] * nm_to_km
                 else:
-                    self.gdf.rmw[it] = np.NaN    
+                    self.gdf.rmw[it] = np.nan    
                 # R35
                 if self.gdf.r35_ne[it] > 0.0:
                     self.gdf.r35_ne[it] = self.gdf.r35_ne[it] * nm_to_km
                 else:
-                    self.gdf.r35_ne[it] = np.NaN
+                    self.gdf.r35_ne[it] = np.nan
 
                 if self.gdf.r35_se[it] > 0.0:
                     self.gdf.r35_se[it] = self.gdf.r35_se[it] * nm_to_km
                 else:
-                    self.gdf.r35_se[it] = np.NaN
+                    self.gdf.r35_se[it] = np.nan
 
                 if self.gdf.r35_sw[it] > 0.0:
                     self.gdf.r35_sw[it] = self.gdf.r35_sw[it] * nm_to_km
                 else:
-                    self.gdf.r35_sw[it] = np.NaN
+                    self.gdf.r35_sw[it] = np.nan
 
                 if self.gdf.r35_nw[it] > 0.0:
                     self.gdf.r35_nw[it] = self.gdf.r35_nw[it] * nm_to_km
                 else:
-                    self.gdf.r35_nw[it] = np.NaN
+                    self.gdf.r35_nw[it] = np.nan
 
                 # r50
                 if self.gdf.r50_ne[it] > 0.0:
                     self.gdf.r50_ne[it] = self.gdf.r50_ne[it] * nm_to_km
                 else:
-                    self.gdf.r50_ne[it] = np.NaN
+                    self.gdf.r50_ne[it] = np.nan
 
                 if self.gdf.r50_se[it] > 0.0:
                     self.gdf.r50_se[it] = self.gdf.r50_se[it] * nm_to_km
                 else:
-                    self.gdf.r50_se[it] = np.NaN
+                    self.gdf.r50_se[it] = np.nan
 
                 if self.gdf.r50_sw[it] > 0.0:
                     self.gdf.r50_sw[it] = self.gdf.r50_sw[it] * nm_to_km
                 else:
-                    self.gdf.r50_sw[it] = np.NaN
+                    self.gdf.r50_sw[it] = np.nan
 
                 if self.gdf.r50_nw[it] > 0.0:
                     self.gdf.r50_nw[it] = self.gdf.r50_nw[it] * nm_to_km
                 else:
-                    self.gdf.r50_nw[it] = np.NaN
+                    self.gdf.r50_nw[it] = np.nan
 
                 # r65
                 if self.gdf.r65_ne[it] > 0.0:
                     self.gdf.r65_ne[it] = self.gdf.r65_ne[it] * nm_to_km
                 else:
-                    self.gdf.r65_ne[it] = np.NaN
+                    self.gdf.r65_ne[it] = np.nan
 
                 if self.gdf.r65_se[it] > 0.0:
                     self.gdf.r65_se[it] = self.gdf.r65_se[it] * nm_to_km
                 else:
-                    self.gdf.r65_se[it] = np.NaN
+                    self.gdf.r65_se[it] = np.nan
 
                 if self.gdf.r65_sw[it] > 0.0:
                     self.gdf.r65_sw[it] = self.gdf.r65_sw[it] * nm_to_km
                 else:
-                    self.gdf.r65_sw[it] = np.NaN
+                    self.gdf.r65_sw[it] = np.nan
 
                 if self.gdf.r65_nw[it] > 0.0:
                     self.gdf.r65_nw[it] = self.gdf.r65_nw[it] * nm_to_km
                 else:
-                    self.gdf.r65_nw[it] = np.NaN
+                    self.gdf.r65_nw[it] = np.nan
 
                 # r100
                 if self.gdf.r100_ne[it] > 0.0:
                     self.gdf.r100_ne[it] = self.gdf.r100_ne[it] * nm_to_km
                 else:
-                    self.gdf.r100_ne[it] = np.NaN
+                    self.gdf.r100_ne[it] = np.nan
 
                 if self.gdf.r100_se[it] > 0.0:
                     self.gdf.r100_se[it] = self.gdf.r100_se[it] * nm_to_km
                 else:
-                    self.gdf.r100_se[it] = np.NaN
+                    self.gdf.r100_se[it] = np.nan
 
                 if self.gdf.r100_sw[it] > 0.0:
                     self.gdf.r100_sw[it] = self.gdf.r100_sw[it] * nm_to_km
                 else:
-                    self.gdf.r100_sw[it] = np.NaN
+                    self.gdf.r100_sw[it] = np.nan
 
                 if self.gdf.r100_nw[it] > 0.0:
                     self.gdf.r100_nw[it] = self.gdf.r100_nw[it] * nm_to_km
                 else:
-                    self.gdf.r100_nw[it] = np.NaN
+                    self.gdf.r100_nw[it] = np.nan
 
             # Done, so set variable
             self.unit_intensity = "ms"
@@ -533,7 +533,7 @@ def add_point(self, time: datetime,
             point = Point(lon, lat)
             gdf_point["geometry"] = [point]
 
-        gdf_point.replace(-999.0, np.NaN, inplace=True)
+        gdf_point.replace(-999.0, np.nan, inplace=True)
         crs = self.gdf.crs
 
         if index == 0:
@@ -970,25 +970,25 @@ def read_trk(filename):
                 gdf_point = GeoDataFrame()
                 gdf_point["geometry"] = [Point(x, y)]  # Assign the new geometry
                 gdf_point["datetime"] = newtime.astype("O").strftime("%Y%m%d %H%M%S")
-                gdf_point["vmax"]   = np.NaN
-                gdf_point["pc"]     = np.NaN
-                gdf_point["rmw"]    = np.NaN
-                gdf_point["r35_ne"] = np.NaN
-                gdf_point["r35_se"] = np.NaN
-                gdf_point["r35_sw"] = np.NaN
-                gdf_point["r35_nw"] = np.NaN
-                gdf_point["r50_ne"] = np.NaN
-                gdf_point["r50_se"] = np.NaN
-                gdf_point["r50_sw"] = np.NaN
-                gdf_point["r50_nw"] = np.NaN
-                gdf_point["r65_ne"] = np.NaN
-                gdf_point["r65_se"] = np.NaN
-                gdf_point["r65_sw"] = np.NaN
-                gdf_point["r65_nw"] = np.NaN
-                gdf_point["r100_ne"] = np.NaN
-                gdf_point["r100_se"] = np.NaN
-                gdf_point["r100_sw"] = np.NaN
-                gdf_point["r100_nw"] = np.NaN
+                gdf_point["vmax"]   = np.nan
+                gdf_point["pc"]     = np.nan
+                gdf_point["rmw"]    = np.nan
+                gdf_point["r35_ne"] = np.nan
+                gdf_point["r35_se"] = np.nan
+                gdf_point["r35_sw"] = np.nan
+                gdf_point["r35_nw"] = np.nan
+                gdf_point["r50_ne"] = np.nan
+                gdf_point["r50_se"] = np.nan
+                gdf_point["r50_sw"] = np.nan
+                gdf_point["r50_nw"] = np.nan
+                gdf_point["r65_ne"] = np.nan
+                gdf_point["r65_se"] = np.nan
+                gdf_point["r65_sw"] = np.nan
+                gdf_point["r65_nw"] = np.nan
+                gdf_point["r100_ne"] = np.nan
+                gdf_point["r100_se"] = np.nan
+                gdf_point["r100_sw"] = np.nan
+                gdf_point["r100_nw"] = np.nan
             else:
                 new_point = False