diff --git a/climaf/plot_operators.py b/climaf/plot_operators.py
index ee551c1..dabbcd5 100644
--- a/climaf/plot_operators.py
+++ b/climaf/plot_operators.py
@@ -32,7 +32,9 @@ def load_plot_operators():
             #
             "--colored_map_file='${in_1}' "
             "--colored_map_variable='${var_1}' "
+            "--colored_map_levels='${clrl}' "
             "--colored_map_levels='${colored_map_levels}' "
+            "--colored_map_cmap='${clrc}' "
             "--colored_map_cmap='${colored_map_color}' "
             "--colored_map_transform=${clrt} "
             "--colored_map_transform_options='${clrto}' "
diff --git a/examples/plotmap_demo.py b/examples/plotmap_demo.py
index 227d013..0cf3c3e 100644
--- a/examples/plotmap_demo.py
+++ b/examples/plotmap_demo.py
@@ -1,12 +1,13 @@
+# # Demonstration of the new CLiMAF operator 'plotmap', dedicated to replace operator 'plot' , which backend script gplot.ncl is at threat due to Ncl maintenance announcements
+
+# ## plotmap is backed by script plotmap.py which uses matplotlib, cartopy and geocat
+
 from IPython.display import Image, display
 from climaf.api import *
-#lcmn
-#export PROJ_DATA=/net/nfs/tools/Users/SU/jservon/spirit-2021.11_envs/20230904/share/proj
-#export PYPROJ_GLOBAL_CONTEXT = ON
-#~/climaf_installs/climaf_running/bin/climaf
 
 
 # +
+# Convenience functions
 def currently_running_in_a_notebook():
     try:
         name = get_ipython().__class__.__name__
@@ -24,14 +25,18 @@ def cshow(obj, drop=False):
         display(Image(cfile(obj)))
     else:
         climaf.api.cshow(obj)
-            
 
 
 # -
 
+# Some CMIP data
 data = fds(
     cpath + "/../examples/data/tas_Amon_CNRM-CM5_historical_r1i1p1_185001-185212.nc")
 
+# More example data
+dg = ds(project="example", simulation="AMIPV6ALB2G",
+        variable="tas", period="1981", frequency="monthly")
+
 # ## Basics
 
 # For reference, a figure using the old script (gplot.ncl)
@@ -40,14 +45,14 @@ def cshow(obj, drop=False):
 # Same with new script. Mapping values to colors is slightly different 
 cshow(plotmap(data, title="Using plotmap"),True)
 
+# One can easily tune the number of levels used
+cshow(plotmap(data, title="hand-tuned levels #", levels=10),True)
+
 # Changing the Reference longitude
 cshow(plotmap(data, title="Changing the Reference longitude", 
               proj="PlateCarree", proj_options={"central_longitude":90},
              ),True)
 
-# One can easily tune the number of levels used
-cshow(plotmap(data, title="hand-tuned levels #", levels=10),True)
-
 # The part of the coordinates range shown is called the 'extent'
 # You can change it, but you must specify longitudes that fits in central_longitude +/- 180°
 # So, by default (i.e. central_longitude=180), that fits in  [ 0, 360 [
@@ -56,36 +61,36 @@ def cshow(obj, drop=False):
 # ## Playing with colors
 
 # Play with number of colors
-cshow(plotmap(data, clrl=15))
+cshow(plotmap(data, levels=15))
 
+# Choosing each color
 mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue,lightseagreen,PaleGreen,Wheat,Brown,Pink'
 cshow(plotmap(data, color= mycolors),True)
 
+# Choosing the mapping of values to colors
 mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue'
 mylevels = '210,220,240,300,310'
 cshow(plotmap(data, color= mycolors, levels=mylevels),True)
 
-mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue'
+# 'colors' and levels" can also be python lists
+mycolors = [ 'white' , 'black' , 'RoyalBlue' , 'LightSkyBlue' , 'PowderBlue' ]
 mylevels = [210,220,240,300,310]
 cshow(plotmap(data, color= mycolors, levels=mylevels),True)
 
-clog('warning')
-mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue'
-mylevels = [210,220,240,300,310]
-cshow(plotmap(data, color= mycolors, levels=mylevels, contours=1),True)
-
-cshow(plotmap(data, color= mycolors, colors=mylevels),True)
+# Adding contours between colored regions
+mylevels = [210,220,240,290,300,310]
+cshow(plotmap(data, levels=mylevels, contours=1),True)
 
+# NCL colormaps are available
 cshow(plotmap(data, color = 'helix'),True)
 
-# Using Ncl colormaps
-cshow(plotmap(nd, clrmap="helix"))
-
+# How to change gridlines look
 cshow(plotmap(data, axis_methods={
     'gridlines': { 'linewidths': 5, 'linestyles': 'dashed', 'draw_labels' : {"bottom": "x", "left": "y"}}}),
       True)
 
-cshow(plotmap(data, clre="contourf"),True)
+# To avoid smoothing, use colormap engine 'pcolormesh' 
+cshow(plotmap(data, clre="pcolormesh"),True)
 
 # ## Data selection
 
@@ -98,23 +103,24 @@ def cshow(obj, drop=False):
 # Selection on date
 cshow(plotmap(data, date='185003', print_time=True),True)
 
-# Selection by providing explicit 'selection options' argument
-cshow(plotmap(data, clrso=[["isel", {"time": 0}]]))
+# Selection by providing explicit 'selection options' arguments
+cshow(plotmap(data, clrso={"isel" : {"time": 0}}, print_time=True),True)
 
 # Again: selection by providing explicit 'selection options' argument
-cshow(plotmap(data, clrso=[['sel', {'time': ' 1850-01'}]]))
+cshow(plotmap(data, clrso={'sel' : {'time': '1850-02'}}, print_time=True))
 
-# ## Mimicking gplot
+# ## Mimicking various gplot options/args
 
 # Changing image size
-cshow(plotmap(data, title="Using plotmap", resolution="600x400"),True)
+cshow(plotmap(data, title="Changing image size", resolution="600x400"),True)
 
 # Showing only a data range. Here is gplot.ncl reference
-cshow(plot(data, min=210, max=250, delta=5),True)
+cshow(plot(data, min=210, max=250, delta=5, title="A range, with 'plot'"),True)
 
 # Showing only a data range. Here is potmap version
-cshow(plotmap(data, min = 210, max = 250, delta=5),True)
+cshow(plotmap(data, min = 210, max = 250, delta=5, title="A range, with 'plotmap'"),True)
 
+# What if surrounding white space is left
 cshow(plotmap(data,trim=False),True)
 
 # Request an horizontal colorbar
@@ -127,15 +133,88 @@ def cshow(obj, drop=False):
 cshow(plotmap(data,focus='land'),True)
 
 # Apply offset and change displayed units accordingly
-cshow(plotmap(data,offset=-273.15, units="C"),True)
+cshow(plotmap(data, offset=-273.15, units="C"),True)
 
 # Polar stereo projection (yet with square limits)
 cshow(plotmap(data,proj='NH70'),True)
 
 # Polylines
-cshow(plotmap(data,xpolyline="0 90",ypolyline="0 45", polyline_options={'color': 'red'}),True)
+cshow(plotmap(data,xpolyline="0 90",ypolyline="0 45", polyline_options={'color': 'blue'}),True)
 
-# ## Playing with projections
+# Nemo data :
+nd = fds(cpath + "/../examples/data/tos_Omon_CNRM_gn_185001-185003.nc")
+cshow(plotmap(nd),True)
+
+# Filling the gaps near ORCA grid poles needs another plot engine  :
+cshow(plotmap(nd, clre="pcolormesh"),True)
+
+# ## Shading
+
+# Representing binary fields using shading. 
+# Binary field #1 is 5th plotmap arg, #2 is 6th arg
+# shdh and shdh2 provides patterns for the binary values
+dgsup=ccdo(dg, operator = "gec,282")
+dginf=ccdo(dg, operator = "lec,252")
+cshow(plotmap(dg, "", "", "", dgsup, dginf,
+              shdh = [None, '/'], shd2h = [None, '+']),True)
+
+# Shading, more dense
+cshow(plotmap(dg, "", "", "", dgsup, dginf,
+      shdh = [None, '///'], shd2h = [None, '+++']),True)
+
+# ## Vectors
+
+# Vector plots. Get data
+ua = fds(cpath + "/../examples/data/uas_CNRM-CM6_sample.nc")
+va = fds(cpath + "/../examples/data/vas_CNRM-CM6_sample.nc")
+
+# Default vector representation type is 'quiver' (arrows)
+# Vector components are args 3 and 4 of plotmap
+cshow(plotmap("","",ua,va),True)
+
+# ### Tuning vector grid size and arrows attributes. See [quiver doc](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html#matplotlib.pyplot.quiver) for more options (but don't play with args X,Y, u, v)
+
+cshow(plotmap("","",ua,va,
+              vecg=50, 
+              veco={'color':'blue', 'headwidth':2.5, 'headlength':4}),
+      True)
+
+# Controling the map projection
+cshow(plotmap("","",ua,va, proj="Mollweide"),True)
+
+# ###  Another vector representation type is barbs. See [its doc](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.barbs.html#matplotlib.pyplot.barbs) for more options
+
+cshow(plotmap("","",ua,va,vecty="barbs", vecg=40, veco={'length':3.5, 'barbcolor' : 'blue'},),True)
+
+# ### The third vector representation type is streamlines. Here is the [doc](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.streamplot.html#matplotlib.pyplot.streamplot) for more options
+
+cshow(plotmap("","",ua,va,vecty="streamplot", veco={'density' : 3 , 'linewidth': 0.7},),True)
+
+# Superimposing colored map and streamplot
+cshow(plotmap(dg,"",ua,va,
+              vecty="streamplot",  
+              veco={'density' : 3 , 'linewidth': 0.7},
+              axis_methods={'set_extent': {'extents': (0. , 60, 0., 60.) }}
+             ),True)
+
+# ## Output formats
+
+# Output in pdf format
+f=cfile(plotmap(dg, format='pdf'))
+print(f)
+# ! display $f
+
+# Output in eps format
+f=cfile(plotmap(dg, format='eps'))
+print(f)
+# ! display $f
+
+# Launching plotmap with format='show' uses 
+# matplotlib.show(block=True), which pops up a window. 
+# This works only outisde a Notebook 
+cfile(plotmap(dg, format="show"))  
+
+# ## Playing with projections, and using data on a rectangular projected grid
 
 # Next data was generated on a grid that uses Lambert2 projection
 l2 = fds(cpath + "/../examples/data/sfcWind_aladin_ext.nc")
@@ -145,9 +224,7 @@ def cshow(obj, drop=False):
 cshow(plotmap(l2),True)
 
 
-# Changing the colored map engine
-cshow(plotmap(l2,clre="pcolormesh"),True)
-
+# This plot engine option speeds up computaton, but may damage the plot
 cshow(plotmap(l2, clreo={'transform_first':True}),True)
 
 # We can request another target projection
@@ -172,7 +249,7 @@ def cshow(obj, drop=False):
 
 
 # If we want to plot on the data grid, we have to describe it explicitly, 
-# using args 'proj' and 'proj_options'. We pack them in a small dict
+# Here we use args 'proj' and 'proj_options', and pack them in a small dict
 projection = { 
     'proj' : 'LambertConformal', 
     'proj_options' : {
@@ -184,44 +261,23 @@ def cshow(obj, drop=False):
 # Note : data will actually be remapped on its own grid (To be checked)
 cshow(plotmap(l2, **projection),True)
 
-# We can also describe target projetion using a file
+# We can also describe target projection using a file with relevant metadata
 cshow(plotmap(l2, proj=cfile(l2)),True)
 
-# Same without data remapping.
+# We can explictly request that the data is not remapped. But there is no watch dog here !
 cshow(plotmap(l2, proj=cfile(l2), clrt="do not remap"),True)
 
-# 
-cshow(plotmap(l2, axis_methods={'set_extent': {'extents': (-60 , 90, -20, 70) }}),True)
-
-# 
-cshow(plotmap(l2, axis_methods={'set_extent': {'extents': (-20 , 80, -20, 70) }},
-             proj="PlateCarree", proj_options={'central_longitude' : 180}),True)
+# Still without remapping, and using another plot engine to 'see' data grid cells
+cshow(plotmap(l2, proj=cfile(l2), clrt="do not remap", clre="pcolormesh"),True)
 
-# time(cshow(plotmap(l2, proj="LambertII", proj_options=[ 10, 37, 37], clre="contourf")))
-# Wall time: 8.82 s
-# time(cshow(plotmap(l2, proj="LambertII", proj_options=[ 10, 37, 37], clre="contourf", clrt="do not remap")))
-# Wall time: 7.48 s
+# Selecting the plot domain, using 'extents'
+cshow(plotmap(l2, axis_methods={'set_extent': {'extents': (0 , 70, 20, 60) }}),True)
 
+# ## Advanced features
 
-# Nemo data :
-nd = fds(cpath + "/../examples/data/tos_Omon_CNRM_gn_185001-185003.nc")
-cshow(plotmap(nd),True)
-
-# Filling the gaps near ORCA grid poles needs another plot engine  :
-cshow(plotmap(nd, clre="pcolormesh"),True)
-
-# Example data
-dg = ds(project="example", simulation="AMIPV6ALB2G",
-        variable="tas", period="1981", frequency="monthly")
-
-cfile(plotmap(dg))
-cfile(plotmap(dg, format="show"))  # Use plt.show(block=True)
-
+# ### The methods of [cartopy.mpl.geoaxes.GeoAxes](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy-mpl-geoaxes-geoaxes) can be called using arg axis_methods
 #
-cfile(plotmap(dg, format='pdf'))  # cshow problématique
-cfile(plotmap(dg, format='eps'))  # cshow problématique
 
-# Axis_methods
 cshow(plotmap(dg, axis_methods={
     'coastlines': {'color': 'grey'},
     'gridlines': {}
@@ -231,164 +287,9 @@ def cshow(obj, drop=False):
 cshow(plotmap(dg, axis_methods={
     'add_feature': {'feature': 'LAND', 'facecolor': 'black', 'zorder': 1}}))
 
-# Pyplot_methods
-# Title, called only once
-cshow(plotmap(dg, plt_methods={
-    'title': {'label': "mytitle", 'loc': 'right'}}))
-
-# Title, with a list of calls
-cshow(plotmap(dg, plt_methods={
-    'title': [
-        {'label': "mytitle", 'loc': 'right'},
-        {'label': "my left title", 'loc': 'left'}
-    ]
-}))
-
-# For pyplot methods with non-keyword args, such as 'plot', use keyword 'largs'
-cshow(plotmap(data, plt_methods={
-    'plot': [{'largs': [[0, 90], [0, 45]], 'color': 'blue', 'marker': 'o'}]
-}),True)
-
+# ### Some methods of [matplotlib.pyplot](https://matplotlib.org/stable/api/pyplot_summary.html#module-matplotlib.pyplot) can be called using arg plt_methods
+#
 
 cshow(plotmap(dg, plt_methods={
     'text': {'x':-120, 'y': 45, 's': 'mytext', 'horizontalalignment': 'left'}}))
 
-
-# ## Shading
-
-# Shading
-dgsup=ccdo(dg, operator = "gec,282")
-dginf=ccdo(dg, operator = "lec,252")
-cshow(plotmap(dg, "", "", "", dgsup, dginf,
-      shdh = [None, '/'], shd2h = [None, '+']),True)
-
-# Shading, more dense
-cshow(plotmap(dg, "", "", "", dgsup, dginf,
-      shdh = [None, '///'], shd2h = [None, '+++']),True)
-
-# ## Vectors
-
-# Vector plots. Get data
-ua = fds(cpath + "/../examples/data/uas_CNRM-CM6_sample.nc")
-va = fds(cpath + "/../examples/data/vas_CNRM-CM6_sample.nc")
-
-cshow(plotmap("","",ua,va),True)
-
-# Default vector type is 'quiver'
-cshow(plotmap("","",ua,va),True)
-
-# Tuning vector grid size and arrows attributes
-cshow(plotmap("","",ua,va,
-              vecg=50, 
-              veco={'color':'blue', 'headwidth':2.5, 'headlength':4}),
-      True)
-
-cshow(plotmap("","",ua,va, proj="Mollweide", vect="PlateCarree"),True)
-
-# +
-
-cshow(plotmap("","",ua,va,vecty="barbs", vecg=40, veco={'length':3.5, 'barbcolor' : 'blue'},),True)
-
-# +
-
-cshow(plotmap(data,"",ua,va,vecty="barbs", vecg=40, veco={'length':3.5, 'barbcolor' : 'blue'},),True)
-# -
-
-# Ploting stremalines yet has a bug re. western longitudes
-cshow(plotmap("","",ua,va,vecty="streamplot", veco={'density' : 3 , 'linewidth': 0.7},),True)
-
-# Ploting stremalines yet has a bug re. western longitudes
-cshow(plotmap("","",ua,va,
-              vecty="streamplot",  
-              veco={'density' : 3 , 'linewidth': 0.7},
-              axis_methods={'set_extent': {'extents': (-180. , 180, -70., 70.) }}
-             ),True)
-
-# +
-
-cshow(plotmap(data,"",ua,va,vecty="streamplot",  veco={'density' : 3 , 'linewidth': 0.7}),True)
-# -
-
-# ## Development
-
-import xarray as xr
-
-f=xr.open_dataset(cpath + "/../examples/data/sfcWind_aladin_ext.nc")
-
-import xarray as xr
-
-f=xr.open_dataset(cpath + "/../examples/data/tas_Amon_CNRM-CM5_historical_r1i1p1_185001-185212.nc")
-
-f.tas.__getattr__('long_name')
-
-dir(f.tas)
-
-f.tas.__getattribute__('isel')
-
-if 'time' in f.tas.isel(time=0).coords:
-    print(f.tas.isel(time=0).time.values)
-
-t=f.tas.isel(time=0).time
-
-t.values
-
-a=plot(ccdo(ds('example|AMIPV6ALB2G|ta|198001|global|monthly'),operator='zonmean'),
-     title='1 field cross-section (without contours lines)',
-     xpolyline='-60.0, -30.0, -30.0, -60.0, -60.0',
-     y='log',
-     ypolyline='70.0, 70.0, 50.0, 50.0, 70.0')
-cshow(a,True)
-
-a=plot(ds('example|AMIPV6ALB2G|tas|1980|global|monthly'),
-       llbox(ds('example|AMIPV6ALB2G|tas|1980|global|monthly'),
-             latmax=80,latmin=30,lonmax=120,lonmin=60),
-       ds('example|AMIPV6ALB2G|uas|1980|global|monthly'),
-       ds('example|AMIPV6ALB2G|vas|1980|global|monthly'),
-       date=19800131,level=10.0,title='Selecting level and time close to 10 and 19800131 respectively',
-       vcRefLengthF=0.02,vcRefMagnitudeF=11.5)
-cshow(a,True)
-
-
-
-cshow(curves(space_average(ds('example|AMIPV6ALB2G|tas|1980-1981|global|monthly')),title='AMIPV6'), True)
-
-from IPython.display import Image, display
-from climaf.api import *
-#lcmn
-#export PROJ_DATA=/net/nfs/tools/Users/SU/jservon/spirit-2021.11_envs/20230904/share/proj
-#export PYPROJ_GLOBAL_CONTEXT = ON
-#~/climaf_installs/climaf_running/bin/climaf
-
-
-# +
-def currently_running_in_a_notebook():
-    try:
-        name = get_ipython().__class__.__name__
-        if name == 'ZMQInteractiveShell':
-            return True
-        else:
-            return False
-    except:
-        return False
-    
-def cshow(obj, drop=False):
-    if drop :
-        cdrop(obj)
-    if currently_running_in_a_notebook():
-        display(Image(cfile(obj)))
-    else:
-        climaf.api.cshow(obj)
-            
-
-
-# -
-
-cshow(plot(ccdo(ccdo(ds('example|AMIPV6ALB2G|ta|1980|global|monthly'),
-                     operator='zonmean'),
-                operator='mermean'),
-           ccdo(
-               ccdo(llbox(ds('example|AMIPV6ALB2G|ta|1980|global|monthly'),
-                          latmax=90,latmin=10,lonmax=150,lonmin=50),operator='zonmean'),
-               operator='mermean'),
-           invXY=True,title='Profiles (t,z)',y='log', resolution='300x300', trim=True,
-          color="BlueDarkRed18"), True)
diff --git a/scripts/plotmap.py b/scripts/plotmap.py
index 26ec53a..10ef817 100644
--- a/scripts/plotmap.py
+++ b/scripts/plotmap.py
@@ -34,6 +34,9 @@ def read_dataset(input_file, variable=None):
 
 
 def filter_dataset(dataset, dimensions, selection_options=dict()):
+    if debug:
+        print("Filter_dataset. selection_options=",
+              selection_options, " dimensions=", dimensions)
     # e.g. selection_options ={ "isel": {"time": 0, "level" : 3} }
     for method in selection_options:
         # e.g. kwargs = {"time": 0, "level" : 3}
@@ -79,11 +82,13 @@ def find_data_in_dataset(variable_dataset, dimensions):
 def find_ccrs(crs_name, options=dict(), data_filename=None):
     """
     Returns a cartopy Coordinate Reference System based on its name and options.
-    If CRS_NAME is a filename, first tries to find both CRS name and options 
-    in file's NetCDF metadata; 
+    If CRS_NAME is a filename, first tries to find both CRS name and options
+    in file's NetCDF metadata;
     Else, if CRS_NAME is None, tries the same with file DATA_FILENAME
     If both unsuccessful, return cartopy.crs.PlateCarree()
     """
+    if debug:
+        print("Find_ccrs with ", crs_name, options, data_filename)
     # First try to identify a file using first crs_name then data_filename
     fic = None
     if crs_name is not None and os.path.exists(crs_name):
@@ -113,10 +118,10 @@ def find_ccrs(crs_name, options=dict(), data_filename=None):
 def ccrs_from_metadata(filename):
     """
     Returns a Cartopy's Coordinate Reference System name and options dict
-    by analyzing NetCDF metadata, assuming it follows the CF convention at 
+    by analyzing NetCDF metadata, assuming it follows the CF convention at
     http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#appendix-grid-mappings
 
-    Yet limited to the case of Lambert Conformal projection, as coded in Aladin model outputs. 
+    Yet limited to the case of Lambert Conformal projection, as coded in Aladin model outputs.
     """
     ccrs_name = None
     ccrs_options = {}
@@ -228,8 +233,8 @@ def get_variable_and_coordinates_from_dataset(
             print(f"Using 2D variables {dimensions} for computing X and Ys")
         lon2d = variable_dataset.coords[dimensions[0]]
         lat2d = variable_dataset.coords[dimensions[1]]
-        #lon2d_fixed = fix_longitudes(lon2d)
-        #X, Y = compute_xy(lat2d, lon2d_fixed, projection, exact)
+        # lon2d_fixed = fix_longitudes(lon2d)
+        # X, Y = compute_xy(lat2d, lon2d_fixed, projection, exact)
         X, Y = compute_xy(lat2d, lon2d, projection, regular)
         variable_coordinates_data.append(X)
         variable_coordinates_data.append(Y)
@@ -335,11 +340,11 @@ def plot_colored_map(fig, ax, coordinates, colored_map_file, colored_map_variabl
     # Apply scale and offset
     if colored_map_scale:
         variable_data = variable_data * colored_map_scale
-        if not arg.units:
+        if not args.units:
             units = "?"
     if colored_map_offset:
         variable_data = variable_data + colored_map_offset
-        if not arg.units:
+        if not args.units:
             units = "?"
 
     # Allow cmap name to refer to package 'cmaps' - see https://github.com/hhuangwx/cmaps
@@ -365,40 +370,42 @@ def plot_colored_map(fig, ax, coordinates, colored_map_file, colored_map_variabl
     contourf_args.update(colored_map_engine_options)
     if colored_map_levels is not None and 'levels' not in contourf_args:
         contourf_args['levels'] = colored_map_levels
-
+    contourf_args["transform"] = transform
+    #
+    if colored_map_engine == "contourf":
+        ploter = plt.contourf
+    elif colored_map_engine == 'pcolormesh':
+        ploter = plt.pcolormesh
+    else:
+        raise ValueError("Unknown plot engine %s" % colored_map_engine)
+    #
     if colored_map_transform != "do not remap":
-        contourf_args["transform"] = transform
-        if colored_map_engine == "contourf":
-            if debug:
-                print("Plotting with contourf and remapping")
-                print("Contourf levels=", colored_map_levels)
-                print("Contourf args=", contourf_args)
-            colored_map_plot = plt.contourf(
-                *variable_coordinates_data, variable_data,
-                **contourf_args)  # colors=colored_map_cmap,
-        elif colored_map_engine == 'pcolormesh':
-            if debug:
-                print("Plotting with pcolormesh and remapping")
-            # if colored_map_cmap is not None:
-            #    raise ValueError("Cannot apply desired colors to the color_map plot "
-            #                     "because plot engine is not set to contourf")
-            norm = create_norm(colored_map_levels, colored_map_cmap,
-                               variable_data, colored_map_min, colored_map_max)
-            colored_map_plot = ax.pcolormesh(
-                *variable_coordinates_data, variable_data, norm=norm, **contourf_args)
-        else:
-            raise ValueError("Unknown plot engine %s" % colored_map_engine)
+        if debug:
+            print(f"Plotting with {colored_map_engine} and remapping; kw=",
+                  contourf_args)
+        if colored_map_engine == 'pcolormesh':
+            contourf_args['norm'] = create_norm(colored_map_levels,
+                                                colored_map_cmap, variable_data,
+                                                colored_map_min, colored_map_max)
+        colored_map_plot = ploter(
+            *variable_coordinates_data, variable_data, **contourf_args)
     else:
         if debug:
-            print("Plotting without remapping, and with contourf")
-        ymin = np.min(variable_coordinates_data[1])
-        ymax = np.max(variable_coordinates_data[1])
-        xmin = np.min(variable_coordinates_data[0])
-        xmax = np.max(variable_coordinates_data[0])
-        contourf_args["extent"] = (xmin, xmax, ymin, ymax)
-        contourf_args["transform"] = projection
-        # Provided colored_map_transform calls for avoiding to remap data
-        colored_map_plot = plt.contourf(variable_data, **contourf_args)
+            print(f"Plotting with {colored_map_engine} and no remapping; kw=",
+                  contourf_args)
+        if colored_map_engine == 'contourf':
+            ymin = np.min(variable_coordinates_data[1])
+            ymax = np.max(variable_coordinates_data[1])
+            xmin = np.min(variable_coordinates_data[0])
+            xmax = np.max(variable_coordinates_data[0])
+            contourf_args["extent"] = (xmin, xmax, ymin, ymax)
+            #contourf_args["transform"] = projection
+            colored_map_plot = plt.contourf(variable_data, **contourf_args)
+        elif colored_map_engine == 'pcolormesh':
+            colored_map_plot = plt.pcolormesh(
+                *variable_coordinates_data, variable_data, **contourf_args)
+
+    #
     for method in colored_map_methods:
         colored_map_plot.__getattr__(method).__call__(
             **colored_map_methods[method])
@@ -434,7 +441,7 @@ def plot_contours_map(ax, coordinates, contours_map_file, contours_map_variable,
                       contours_map_selection_options, contours_map_min, contours_map_max,
                       contours_map_scale, contours_map_offset, title, title_is_plotted):
     # Find the transform
-    if contours_map_transform:
+    if contours_map_transform != "do not remap":
         transform = find_ccrs(contours_map_transform,
                               contours_map_transform_options,
                               contours_map_file)
@@ -470,11 +477,11 @@ def plot_contours_map(ax, coordinates, contours_map_file, contours_map_variable,
         kwargs["vmax"] = contours_map_max
     if len(contours_map_colors) > 0:
         kwargs["colors"] = contours_map_colors
-    if contours_map_levels is not None:
-        plt.contour(*variable_coordinates_data, variable_data,
-                    contours_map_levels, **kwargs)
-    else:
-        plt.contour(*variable_coordinates_data, variable_data, **kwargs)
+    if contours_map_levels is not None and 'levels' not in kwargs:
+        kwargs['levels'] = contours_map_levels
+    if debug:
+        print("Plotting with contour and args", kwargs)
+    plt.contour(*variable_coordinates_data, variable_data, **kwargs)
 
     if not title_is_plotted:
         gv.set_titles_and_labels(
@@ -487,7 +494,7 @@ def plot_shaded_map(ax, coordinates, shaded_map_file, shaded_map_variable,
                     shaded_map_min, shaded_map_max,
                     shaded_map_scale, shaded_map_offset, title, title_is_plotted):
     # Find the transform
-    if shaded_map_transform:
+    if shaded_map_transform != "do not remap":
         transform = find_ccrs(shaded_map_transform,
                               shaded_map_transform_options,
                               shaded_map_file)
@@ -534,7 +541,7 @@ def plot_vector_map(ax, coordinates, vectors_map_u_file, vectors_map_v_file,
                     vectors_map_type, vectors_map_options, vectors_map_selection_options,
                     vectors_map_scale, vectors_map_gridsizes, title, title_is_plotted):
     # Find the transform
-    if vectors_map_transform:
+    if vectors_map_transform != "do not remap":
         transform = find_ccrs(vectors_map_transform,
                               vectors_map_transform_options,
                               vectors_map_u_file)
@@ -610,7 +617,8 @@ def plot_map(args, polar_stereo_extent=None, print_time=False):
     # Deal with projection
     projection = find_ccrs(args.projection, args.projection_options)
     if debug:
-        print("Map projection set to %s" % args.projection)
+        print("Map projection set to %s %s" %
+              (args.projection, args.projection_options))
 
     # Initialize the plot
     fig = plt.figure(**args.figure_options)
@@ -628,7 +636,7 @@ def plot_map(args, polar_stereo_extent=None, print_time=False):
         for kwargs in args.axis_methods[method]:
             largs = kwargs.pop('largs', [])
             if method == 'set_extent':
-                #kwargs['crs'] = projection
+                # kwargs['crs'] = projection
                 kwargs['extents'] = tuple(kwargs['extents'])
             if debug:
                 print(f"Calling axis method {method} with kwargs: {kwargs}")
diff --git a/scripts/plotmap_parsing.py b/scripts/plotmap_parsing.py
index 3b047e3..fa2d129 100644
--- a/scripts/plotmap_parsing.py
+++ b/scripts/plotmap_parsing.py
@@ -374,7 +374,7 @@ def process_args(args):
                         if feature not in features:
                             raise ValueError(f"Required feature {feature} doesn't " +
                                              f"belong to cartopy.features {features}")
-                        dic['feature'] = eval(f"cartopy.cfeature.{feature}")
+                        dic['feature'] = eval(f"cartopy.feature.{feature}")
                 else:
                     raise ValueError(
                         "axis_methods 'add_feature' does not include key 'feature' in its args  dic :",
@@ -451,6 +451,32 @@ def mimic_gplot(args, selection_options_list):
         args.savefig_options['bbox_inches'] = 'tight'
         # else just let what the caller may have set
 
+    # Projection shorcuts
+    if args.projection is not None and args.projection[0:2] in ['NH', 'SH']:
+        print(
+            "TBD : For proj == NHxx or SHxx, limit is not yet a circle. " +
+            "This can be improved using that URL: " +
+            "https://scitools.org.uk/cartopy/docs/latest/gallery/lines_and_polygons/always_circular_stereo.html")
+        if len(args.projection) > 2:
+            latitude_limit = float(args.projection[2:])
+        if args.projection[0:2] == 'NH':
+            args.projection = "NorthPolarStereo"
+            settings['polar_stereo_extent'] = [-180, 180, latitude_limit, 90]
+        if args.projection[0:2] == 'SH':
+            args.projection = "SouthPolarStereo"
+            settings['polar_stereo_extent'] = [-180, 180, -90, -latitude_limit]
+        if args.projection_options is None:
+            args.projection_options = dict()
+        if "central_longitude" not in args.projection_options:
+            args.projection_options["central_longitude"] = 0.0
+
+    if args.projection_options is None:
+        if args.projection is None:
+            args.projection = "PlateCarree"
+            args.projection_options = {'central_longitude': 180.}
+        else:
+            args.projection_options = {}
+
     # Coastlines. Set it by default, and allow user to override default
     if 'coastlines' not in args.axis_methods:
         args.axis_methods['coastlines'] = [{}]
@@ -492,9 +518,12 @@ def mimic_gplot(args, selection_options_list):
     #
     #
     if args.colored_map_cmap is not None:
-        if "," in args.colored_map_cmap:
+        if "," in args.colored_map_cmap or type(args.colored_map_cmap) is list:
             # We have an explicit list of color names -> use contourf args 'colors'
-            cdic['colors'] = args.colored_map_cmap.split(",")
+            if "," in args.colored_map_cmap:
+                cdic['colors'] = args.colored_map_cmap.split(",")
+            else:
+                cdic['colors'] = args.colored_map_cmap
             cdic['cmap'] = None
             # In that case contourf doesn't support vmin/vmax, nor 'norm'
             # and maybe needs 'levels'.
@@ -522,11 +551,11 @@ def mimic_gplot(args, selection_options_list):
     if args.debug:
         print("cdic=", cdic)
     args.colored_map_engine_options.update(**cdic)
-    if args.contours_map_levels == 1 and args.colored_map_levels is not None:
-        args.contours_map_levels = args.colored_map_levels
 
-    # contours = 1 allows to simply draw contours of the colored map field by
-    # providing contours_map_levels (a check is done that no contours file is provided)
+    # arg 'contours_map_level' (or (contours')' allows to simply draw
+    # contours of the colored map field by providing
+    # contours_map_levels (a check is done that no contours file is
+    # provided)
     if args.contours_map_levels is not None and \
        args.contours_map_file is None and \
        args.colored_map_file is not None:
@@ -534,6 +563,8 @@ def mimic_gplot(args, selection_options_list):
         args.contours_map_variable = args.colored_map_variable
         args.contours_map_transform = args.colored_map_transform
         args.contours_map_transform_options = args.colored_map_transform_options
+        if args.contours_map_levels == 1 and args.colored_map_levels is not None:
+            args.contours_map_levels = args.colored_map_levels
         if args.contours_map_colors is None:
             args.contours_map_colors = args.colored_map_cmap
 
@@ -549,7 +580,7 @@ def mimic_gplot(args, selection_options_list):
         if 'add_feature' not in args.axis_methods:
             args.axis_methods['add_feature'] = []
         args.axis_methods['add_feature'].append(
-            {'feature': eval(f"cfeature.{feature}"), 'facecolor': 'white', 'zorder': 1})
+            {'feature': eval(f"cartopy.feature.{feature}"), 'facecolor': 'white', 'zorder': 1})
 
     if args.scale:
         if args.colored_map_file is not None:
@@ -615,32 +646,6 @@ def mimic_gplot(args, selection_options_list):
                 options['sel'] = dict()
             options['sel']['time'] = date
 
-    # Projection shorcuts
-    if args.projection is not None and args.projection[0:2] in ['NH', 'SH']:
-        print(
-            "TBD : For proj == NHxx or SHxx, limit is not yet a circle. " +
-            "This can be improved using that URL: " +
-            "https://scitools.org.uk/cartopy/docs/latest/gallery/lines_and_polygons/always_circular_stereo.html")
-        if len(args.projection) > 2:
-            latitude_limit = float(args.projection[2:])
-        if args.projection[0:2] == 'NH':
-            args.projection = "NorthPolarStereo"
-            settings['polar_stereo_extent'] = [-180, 180, latitude_limit, 90]
-        if args.projection[0:2] == 'SH':
-            args.projection = "SouthPolarStereo"
-            settings['polar_stereo_extent'] = [-180, 180, -90, -latitude_limit]
-        if args.projection_options is None:
-            args.projection_options = dict()
-        if "central_longitude" not in args.projection_options:
-            args.projection_options["central_longitude"] = 0.0
-
-    if args.projection_options is None:
-        if args.projection is None:
-            args.projection = "PlateCarree"
-            args.projection_options = {'central_longitude': 180.}
-        else:
-            args.projection_options = {}
-
     if args.xpolyline is not None and args.ypolyline is not None:
         x = args.xpolyline.split()
         y = args.ypolyline.split()