Psurf

_install/setup_epygram.py

@@ -78,7 +78,7 @@ def vortex_installdir(self):
         install_dirs = {
             'default': '/home/common/sync/vortex/vortex-olive',
             'bullx': '/home/mf/dp/marp/verolive/vortex/vortex-olive',
-            'atos_sequana': '/home/mf/dp/marp/verolive/vortex/vortex-olive',
+            'atos_sequana': '/home/sos/vortex/vortex-olive',
             'dsidev': '/soprano/home/marp999/vortex/vortex-olive',
             'ecmwf_cc': '/home/ms/fr/sos/vortex/vortex-olive',
             'ecgate': '/home/ms/fr/sos/vortex/vortex-olive'}

apptools/epy_conv.py

@@ -17,8 +17,6 @@
 import footprints
 from footprints import FPDict, FPList
 import taylorism
-assert LooseVersion(taylorism.__version__) >= LooseVersion('1.0.7'), \
-    'Version of package taylorism (here {}) need to be >= 1.0.7'.format(taylorism.__version__)
 from bronx.fancies.display import printstatus
 from bronx.syntax.parsing import str2dict
 

docs/source/gallery/C-data/08.morphing_apply_local_deformation.ipynb

@@ -0,0 +1,137 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Morphing: apply a local deformation of field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "# for figures in notebook\n",
+    "\n",
+    "# import & initialize epygram\n",
+    "import epygram\n",
+    "epygram.init_env()\n",
+    "\n",
+    "import os\n",
+    "INPUTS_DIR = os.path.join('..', 'inputs')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r = epygram.formats.resource(os.path.join(INPUTS_DIR, 'ICMSHAROM+0022'), 'r')\n",
+    "f = r.readfield('SURFTEMPERATURE')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = f.cartoplot(subzone='CI')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# make 2 PointFields by extracting them from the field\n",
+    "p1 = f.extract_point(8,45)\n",
+    "p2 = f.extract_point(8.5, 45.5)\n",
+    "print(p1.getdata(), p2.getdata())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set arbitrary different values\n",
+    "p1.setdata(270)\n",
+    "p2.setdata(285)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f.morph_with_points([p1, p2], morphing='gaussian', sigma=10000)  # gaussian blending weight with 10km sigma"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = f.cartoplot(subzone='CI')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# or as increment (less brutal):\n",
+    "p = f.extract_point(7,46)\n",
+    "p.setdata(+10)  # +10K increment\n",
+    "f.morph_with_points([p], morphing='gaussian', sigma=10000, increment=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = f.cartoplot(subzone='CI')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}

docs/source/gallery/index.rst

@@ -54,6 +54,7 @@ C - Data
 #. `Recompute u/v from psi/khi (Arpege) <C-data/05.wind_transformations_psikhi.ipynb>`_
 #. `Compute vor/div from u/v <C-data/06.wind_transformations_vordiv.ipynb>`_
 #. `Shaving: cut values over a threshold down to this threshold <C-data/07.shave.ipynb>`_
+#. `Morphing: apply local deformation <C-data/08.morphing_apply_local_deformation.ipynb>`_
 
 D.1 - Horizontal Geometry
 -------------------------

epygram/__init__.py

@@ -83,10 +83,11 @@
 
 import footprints
 
+package_rootdir = os.path.dirname(os.path.realpath(__path__[0]))  # realpath to resolve symlinks
+
 __all__ = []
 
-__version__ = io.open(os.path.join(os.path.dirname(os.path.dirname(os.path.realpath(__file__))),
-                                   'VERSION'), 'r').read().strip()
+__version__ = io.open(os.path.join(package_rootdir, 'VERSION'), 'r').read().strip()
 
 __license__ = 'CeCILL-C'
 

epygram/_plugins/with_cartopy/H2DField.py

@@ -315,13 +315,21 @@ def _cartoplot_shape(self,
     if plot_method in ('contourf', 'contour'):
         if not self.geometry.rectangular_grid:
             # gauss grid
+            # There are three kinds of masked data:
+            #  - due to the shape of the gaussian grid (already masked in data, x and y):
+            #      => must be suppressed from the output arrays
+            #  - due to the projection (already masked in x and y):
+            #      => must be suppressed from the output arrays
+            #  - due to undefined values (masked in data, eg: SST over continents):
+            #      => must be kept as missing in the output arrays (to actually have a hole on the plot)
             _fill_value = 1e20
-            _shp = [1,1] + list(data.shape)
+            _shp = [1, 1] + list(data.shape)
             zf = self.geometry.fill_maskedvalues(data.reshape(_shp),  # protect actually masked values
-                                                 _fill_value)
+-                                                _fill_value).reshape(data.shape)
+            zf = numpy.ma.masked_where(xf.mask, numpy.ma.masked_where(yf.mask, zf)) # zf masked if xf or yf is masked
+            xf = numpy.ma.array(numpy.ma.masked_where(zf.mask, x).compressed()) # xf masked where zf is masked, then compressed
+            yf = numpy.ma.array(numpy.ma.masked_where(zf.mask, y).compressed()) # yf masked where zf is masked, then compressed
             zf = numpy.ma.masked_equal(zf.compressed(), _fill_value)  # flatten and re-masked
-            xf = numpy.ma.masked_where(zf.mask, x.compressed())
-            yf = numpy.ma.masked_where(zf.mask, y.compressed())
         elif self.geometry.dimensions['Y'] == 1:
             # unstructured grid
             xf = x.flatten()

epygram/fields/D3Field.py

@@ -923,7 +923,7 @@ def extract_zoom(self, zoom, extra_10th=False):
                     zoomvals = []
                     for i in flat_indexes:
                         zoomvals.append(vals[i])
-                    data[t, k, :, :] = numpy.array(zoomvals).reshape(shp[2:])
+                    data[t, k, :, :] = numpy.ma.array(zoomvals).reshape(shp[2:])
 
         fid = {k:v for k, v in self.fid.items()}
         for k, v in fid.items():
@@ -2433,7 +2433,7 @@ def getvalidity(self, index_or_validity):
             # levels are, at least, dependent on position
             levels4d = newfield.geometry.get_levels(d4=True, nb_validities=len(self.validity))
             kwargs_vcoord = copy.deepcopy(newfield.geometry.vcoordinate.footprint_as_dict())
-            kwargs_vcoord['levels'] = levels4d[index, ...].squeeze()
+            kwargs_vcoord['levels'] = list(levels4d[index, ...].squeeze())
             newfield.geometry.vcoordinate = fpx.geometry(**kwargs_vcoord)
 
         return newfield

epygram/fields/H1DField.py

@@ -532,7 +532,7 @@ def update(i, ax, myself, transecti, title_prefix, kwargs):
             ax.clear()
             transecti = myself.getvalidity(i)
             if title_prefix is not None:
-                title = title_prefix + '\n' + transecti.validity.get().isoformat(sep=b' ')
+                title = title_prefix + '\n' + transecti.validity.get().isoformat(sep=' ')
             transecti.plotfield(title=title, **kwargs)
 
     anim = animation.FuncAnimation(fig, update,

epygram/fields/H2DField.py

@@ -104,6 +104,28 @@ def extract_point(self, lon, lat,
         pt.setdata(value)
         return pt
 
+    def extract_contour(self,level):
+        """
+        Return the intersection of the plan **level** and the H2Dfield: the result is a list of (lon,lat) point list.
+        The length of the first list is the number of disjoint lines.
+        This function relies on the contourpy library which is expected to be installed by the user.
+        """
+        try:
+            from contourpy import contour_generator
+        except ModuleNotFoundError:
+            raise ModuleNotFoundError('To use *extract_contour* method, you must have `contourpy` installed: ' +
+                                      'pip3 install --user contourpy')
+        cont_gen=contour_generator(z=self.getdata()) 
+        contour=cont_gen.lines(level)
+        lines=[]
+        for ijlist in contour:
+            line=[]
+            for i,j in ijlist:
+                lon,lat=self.geometry.ij2ll(i,j)
+                line.append((lon,lat))
+            lines.append(line)
+        return lines
+
 ###################
 # PRE-APPLICATIVE #
 ###################

epygram/fields/V1DField.py

@@ -480,7 +480,7 @@ def plotanimation(profile,
             title_prefix = ''
         else:
             title_prefix = title
-        title = title_prefix + '\n' + profile.validity[0].get().isoformat(sep=b' ')
+        title = title_prefix + '\n' + profile.validity[0].get().isoformat(sep=' ')
     else:
         title_prefix = None
     profile0 = profile.getvalidity(0)
@@ -511,7 +511,7 @@ def update(i, ax, myself, profilei, title_prefix, kwargs):
             ax.clear()
             profilei = myself.getvalidity(i)
             if title_prefix is not None:
-                title = title_prefix + '\n' + profilei.validity.get().isoformat(sep=b' ')
+                title = title_prefix + '\n' + profilei.validity.get().isoformat(sep=' ')
             profilei.plotfield(title=title,
                                **kwargs)
 

epygram/fields/V2DField.py

@@ -321,7 +321,7 @@ def plotanimation(self,
                 title_prefix = ''
             else:
                 title_prefix = title
-            title = title_prefix + '\n' + self.validity[0].get().isoformat(sep=b' ')
+            title = title_prefix + '\n' + self.validity[0].get().isoformat(sep=' ')
         else:
             title_prefix = None
         field0 = self.getvalidity(0)
@@ -343,7 +343,7 @@ def update(i, ax, myself, fieldi, title_prefix, kwargs):
                 ax.clear()
                 fieldi = myself.getvalidity(i)
                 if title_prefix is not None:
-                    title = title_prefix + '\n' + fieldi.validity.get().isoformat(sep=b' ')
+                    title = title_prefix + '\n' + fieldi.validity.get().isoformat(sep=' ')
                 fieldi.plotfield(title=title,
                                  **kwargs)
 

epygram/formats/FA.py

@@ -1276,10 +1276,12 @@ def extract_subdomain(self,
                 Psurf = self.readfield('SURFPRESSION')
                 if Psurf.spectral:
                     Psurf.sp2gp()
-                ps_transect = numpy.exp(Psurf.getvalue_ll(*geometry.get_lonlat_grid(),
+                ps_transect = Psurf.getvalue_ll(*geometry.get_lonlat_grid(),
                                                           interpolation=interpolation,
                                                           one=False,
-                                                          external_distance=external_distance))
+                                                          external_distance=external_distance)
+                if ps_transect.mean() < 100. : 
+                    ps_transect = numpy.exp(ps_transect)
                 ps_transect = ps_transect.reshape(h_shape)
                 del Psurf
             # P => H necessary profiles

epygram/geometries/D3Geometry.py

@@ -3033,8 +3033,8 @@ def select_subzone(self, subzone):
                                              'X_Czone':self.dimensions['X_Czone'],
                                              'Y_Czone':self.dimensions['Y_Czone']}
             elif subzone == 'C':
-                centerPoint = ((float(self.dimensions['X']) - 1) / 2.,
-                               (float(self.dimensions['Y']) - 1) / 2.)  # Coordinates of center point
+                centerPoint = ((float(self.dimensions['X_Czone']) - 1) / 2.,
+                               (float(self.dimensions['Y_Czone']) - 1) / 2.)  # Coordinates of center point
                 geom_kwargs['grid']['LAMzone'] = None
                 geom_kwargs['dimensions'] = {'X':self.dimensions['X_Czone'],
                                              'Y':self.dimensions['Y_Czone']}
@@ -4592,8 +4592,8 @@ def _rng(n):
             result = list(result)  # We transform into list to be able to modify the length
             for ipt, points in enumerate(result):
                 if request.get('shape', 'circle') == 'circle':
-                    for i,j in points:
-                        print(i,j,self.distance((lon[ipt], lat[ipt]), self.ij2ll(i, j)))
+                    #for i,j in points:
+                    #    print(i,j,self.distance((lon[ipt], lat[ipt]), self.ij2ll(i, j)))
                     result[ipt] = [(i, j)
                                    for (i, j) in points
                                    if self.distance((lon[ipt], lat[ipt]), self.ij2ll(i, j)) <= request['radius']]

site/epyweb/actions.py

@@ -114,13 +114,12 @@ def POST(self):
             if 'date' in vortexArgs:
                 vortexArgs['date'] = util.datex(vortexArgs['date'].encode())
             if 'term' in vortexArgs:
-                print(vortexArgs['term']) #.encode())
                 vortexArgs['term'] = rangex(vortexArgs['term']) #.encode())
             if 'month' in vortexArgs:
                 vortexArgs['month'] = rangex(vortexArgs['month'].encode())
             if 'member' in vortexArgs:
                 memberstr = '_mb[member]'
-                vortexArgs['member'] = rangex(vortexArgs['member'].encode())
+                vortexArgs['member'] = rangex(vortexArgs['member']) #.encode())
             else:
                 memberstr = ''
 

tests/get_data.sh

@@ -5,7 +5,7 @@ PORT=21
 LOGIN=epygram  # read-only anonymous user
 PASSWORD=EpygraM  # with anonymous password
 
-ftp -i -v -n $HOST $PORT << EOF
+ftp -i -p -v -n $HOST $PORT << EOF
 quote USER $LOGIN
 quote PASS $PASSWORD
 bin