Merge pull request #161 from brey/dev

Dev

pyposeidon/boundary.py

@@ -302,13 +302,17 @@ def tag(geometry, coasts, cbuffer, blevels):
 
             ww = gp.GeoDataFrame(geometry=cs)
 
-            try:
-                gw = gp.GeoDataFrame(
-                    geometry=list(ww.buffer(0).unary_union)
-                )  # merge the polygons that are split (around -180/180)
-            except:
+            if ww.empty:
                 gw = gp.GeoDataFrame(geometry=list(ww.values))
 
+            else:
+                try:
+                    gw = gp.GeoDataFrame(
+                        geometry=list(ww.buffer(0).unary_union)
+                    )  # merge the polygons that are split (around -180/180)
+                except:
+                    gw = gp.GeoDataFrame(geometry=list(ww.values))
+
             if wc.geom_type.all() != "Polygon":
                 gw = gp.GeoDataFrame(geometry=gw.boundary.values)
 

pyposeidon/dem.py

@@ -87,7 +87,7 @@ def __init__(self, dem_source: str, **kwargs):
                 self.adjust(coastline, **kwargs)
 
     def adjust(self, coastline, **kwargs):
-        self.Dataset, check = fix(self.Dataset, coastline, **kwargs)
+        self.Dataset, check, flag = fix(self.Dataset, coastline, **kwargs)
 
         if not check:
             logger.warning("Adjusting dem failed, keeping original values\n")

pyposeidon/schism.py

@@ -566,7 +566,11 @@ def output(self, **kwargs):
                 logger.info("Keeping bathymetry from hgrid.gr3 ..\n")
 
                 copyfile(self.mesh_file, os.path.join(path, "hgrid.gr3"))  # copy original grid file
-                copyfile(os.path.join(path, "hgrid.gr3"), os.path.join(path, "hgrid.ll"))
+                src = os.path.join(path, "hgrid.gr3")
+                dst = os.path.join(path, "hgrid.ll")
+                if os.path.lexists(dst):
+                    os.remove(dst)
+                os.symlink(src, dst)
 
             # manning file
             manfile = os.path.join(path, "manning.gr3")
@@ -677,6 +681,11 @@ def run(self, **kwargs):
 
         proc = tools.execute_schism_mpirun_script(cwd=calc_dir)
 
+        if proc.returncode == 0:
+            # ---------------------------------------------------------------------
+            logger.info("model finished successfully\n")
+            # ---------------------------------------------------------------------
+
     def save(self, **kwargs):
         path = get_value(self, kwargs, "rpath", "./schism/")
 
@@ -705,7 +714,13 @@ def save(self, **kwargs):
                 dic.update({"meteo": [x.attrs for x in meteo.Dataset]})
 
         coastlines = self.__dict__.get("coastlines", None)
-        dic.update({"coastlines": coastlines})
+        if coastlines is not None:
+            # Save the path to the serialized coastlines - #130
+            coastlines_database = os.path.join(path, "coastlines.json")
+            coastlines.to_file(coastlines_database)
+            dic.update({"coastlines": coastlines_database})
+        else:
+            dic.update({"coastlines": None})
 
         dic["version"] = pyposeidon.__version__
 
@@ -1297,10 +1312,16 @@ def results(self, **kwargs):
         path = get_value(self, kwargs, "rpath", "./schism/")
 
         logger.info("get combined 2D netcdf files \n")
+        # check for new IO output
         hfiles = glob.glob(os.path.join(path, "outputs/out2d_*.nc"))
         hfiles.sort()
+        # check for old IO output
+        ofiles = glob.glob(os.path.join(path, "outputs/schout_*_*.nc"))
+
+        if not (ofiles or hfiles):
+            logger.warning("no output netcdf files, moving on")
 
-        if hfiles:
+        elif hfiles:
             x2d = xr.open_mfdataset(hfiles, data_vars="minimal")
 
             # set timestamp
@@ -1335,7 +1356,7 @@ def results(self, **kwargs):
             # save 2D variables to file
             x2d.to_netcdf(os.path.join(path, "outputs/schout_1.nc"))
 
-        else:
+        elif ofiles:
             if len(self.misc) == 0:
                 logger.info("retrieving index references ... \n")
                 self.global2local(**kwargs)
@@ -1704,6 +1725,9 @@ def results(self, **kwargs):
 
                 xc.to_netcdf(os.path.join(path, f"outputs/schout_{val}.nc"))
 
+        else:
+            raise Exception("This should never happen")
+
         logger.info("done with output netCDF files \n")
 
     def set_obs(self, **kwargs):
@@ -1733,6 +1757,9 @@ def set_obs(self, **kwargs):
         ##### normalize to be used inside pyposeidon
         tgn = normalize_column_names(tg.copy())
 
+        ##### make sure lat/lon are floats
+        tgn = tgn.astype({"latitude": float, "longitude": float})
+
         coastal_monitoring = get_value(self, kwargs, "coastal_monitoring", False)
         flags = get_value(self, kwargs, "station_flags", [1] + [0] * 8)
 
@@ -1833,7 +1860,7 @@ def set_obs(self, **kwargs):
         self.params["SCHOUT"]["nspool_sta"] = nspool_sta
         self.params.write(os.path.join(path, "param.nml"), force=True)
 
-        self.stations = stations
+        self.stations_mesh_id = stations
 
         logger.info("write out stations.in file \n")
 
@@ -1853,25 +1880,25 @@ def get_station_sim_data(self, **kwargs):
         try:
             # get the station flags
             flags = pd.read_csv(os.path.join(path, "station.in"), header=None, nrows=1, delim_whitespace=True).T
-            flags.columns = ["flag"]
-            flags["variable"] = [
-                "elev",
-                "air_pressure",
-                "windx",
-                "windy",
-                "T",
-                "S",
-                "u",
-                "v",
-                "w",
-            ]
-
-            vals = flags[flags.values == 1]  # get the active ones
-        except OSError as e:
-            if e.errno == errno.EEXIST:
-                logger.error("no station.in file present")
+        except FileNotFoundError:
+            logger.error("no station.in file present")
             return
 
+        flags.columns = ["flag"]
+        flags["variable"] = [
+            "elev",
+            "air_pressure",
+            "windx",
+            "windy",
+            "T",
+            "S",
+            "u",
+            "v",
+            "w",
+        ]
+
+        vals = flags[flags.values == 1]  # get the active ones
+
         dstamp = kwargs.get("dstamp", self.rdate)
 
         dfs = []

pyposeidon/utils/cast.py

@@ -285,6 +285,10 @@ def run(self, **kwargs):
         self.rdate = self.model.rdate
         ppath = self.ppath
 
+        # ppath = pathlib.Path(ppath).resolve()
+        # ppath = str(ppath)
+        ppath = os.path.realpath(ppath)
+
         # control
         if not isinstance(self.rdate, pd.Timestamp):
             self.rdate = pd.to_datetime(self.rdate)
@@ -296,6 +300,10 @@ def run(self, **kwargs):
         # create the new folder/run path
         rpath = self.cpath
 
+        #    rpath = pathlib.Path(rpath).resolve()
+        #     rpath = str(rpath)
+        rpath = os.path.realpath(rpath)
+
         if not os.path.exists(rpath):
             os.makedirs(rpath)
 
@@ -309,11 +317,16 @@ def run(self, **kwargs):
             info = data.to_dict(orient="records")[0]
 
         try:
-            args = set(kwargs.keys()).intersection(info.keys())  # modify dic with kwargs
+            args = info.keys() & kwargs.keys()  # modify dic with kwargs
             for attr in list(args):
-                info[attr] = kwargs[attr]
-        except:
-            pass
+                if isinstance(info[attr], dict):
+                    info[attr].update(kwargs[attr])
+                else:
+                    info[attr] = kwargs[attr]
+                setattr(self, attr, info[attr])
+        except Exception as e:
+            logger.exception("problem with kwargs integration\n")
+            raise e
 
         # add optional additional kwargs
         for attr in kwargs.keys():
@@ -333,30 +346,14 @@ def run(self, **kwargs):
 
         m = pm.set(**info)
 
-        # Mesh
-        gfile = glob.glob(os.path.join(ppath, "hgrid.gr3"))
-        if gfile:
-            info["mesh_file"] = gfile[0]
-            self.mesh_file = gfile[0]
-            info["mesh_generator"] = None
-            self.mesh_generator = None
-
-        m.mesh = pmesh.set(type="tri2d", **info)
-
-        # get lat/lon from file
-        if hasattr(self, "mesh_file"):
-            info.update({"lon_min": m.mesh.Dataset.SCHISM_hgrid_node_x.values.min()})
-            info.update({"lon_max": m.mesh.Dataset.SCHISM_hgrid_node_x.values.max()})
-            info.update({"lat_min": m.mesh.Dataset.SCHISM_hgrid_node_y.values.min()})
-            info.update({"lat_max": m.mesh.Dataset.SCHISM_hgrid_node_y.values.max()})
-
         # copy/link necessary files
         logger.debug("Copy necessary + station files")
         copy_files(rpath=rpath, ppath=ppath, filenames=self.files + self.station_files)
-        logger.debug("Copy model files")
         if copy:
+            logger.debug("Copy model files")
             copy_files(rpath=rpath, ppath=ppath, filenames=self.model_files)
         else:
+            logger.debug("Symlink model files")
             symlink_files(rpath=rpath, ppath=ppath, filenames=self.model_files)
 
         logger.debug(".. done")
@@ -391,18 +388,13 @@ def run(self, **kwargs):
         else:
             logger.info("Symlinking`: %s -> %s", inresfile, outresfile)
             try:
-                os.symlink(
-                    pathlib.Path(os.path.join(ppath, inresfile)).resolve(strict=True), os.path.join(rpath, outresfile)
-                )
+                os.symlink(inresfile, outresfile)
             except OSError as e:
                 if e.errno == errno.EEXIST:
                     logger.warning("Restart link present\n")
                     logger.warning("overwriting\n")
-                    os.remove(os.path.join(rpath, outresfile))
-                    os.symlink(
-                        pathlib.Path(os.path.join(ppath, inresfile)).resolve(strict=True),
-                        os.path.join(rpath, outresfile),
-                    )
+                    os.remove(outresfile)
+                    os.symlink(inresfile, outresfile)
                 else:
                     raise e
         # get new meteo

pyposeidon/utils/data.py

@@ -192,11 +192,15 @@ def __init__(self, **kwargs):
 
                 datai.append(xdat)  # append to list
 
-        merge = kwargs.get("merge", True)
+        if not any(datai):
+            logger.warning("no output netcdf files.")
+            self.Dataset = None
+        else:
+            merge = kwargs.get("merge", True)
 
-        if merge:
-            datai = flat_list(datai)
-            self.Dataset = xr.open_mfdataset(datai, combine="by_coords", data_vars="minimal")
+            if merge:
+                datai = flat_list(datai)
+                self.Dataset = xr.open_mfdataset(datai, combine="by_coords", data_vars="minimal")
 
-        else:
-            self.Dataset = [xr.open_mfdataset(x, combine="by_coords", data_vars="minimal") for x in datai]
+            else:
+                self.Dataset = [xr.open_mfdataset(x, combine="by_coords", data_vars="minimal") for x in datai]