Re: [ferret_users] restrict polymark to transect

[Marco van Hulten <marco.van-hulten@xxxxxxxxxxxx>]

Le 2015-10-08 Marco a écrit:
> Hello,
> 
> I regularly plot transects of model output with observations plotted
> on top of these as dots.  So the last command is
> 
>   go polymark polygon/over Latitude, Depth, myvar circle 1.0
> 
> for a transect that is monotonic in latitude.  Observational data from
> all longitudes are plotted.  This works fine when using only the
> transect observational data.  If I use a dataset that includes
> observations at longitudes far away from the transect (but within the
> latitude range of the transect), these are also included on the plot.
> 
> While this is expected behaviour, I would like to limit the longitudes
> to close to the transect.  So I am thinking of including a criterium
> that limits plotting of points near the transect.  How would such a
> criterium look like?
> [...]

Paulo and Ryo suggested that I could define a mask for this.  I am
trying to extend this to an arbitrary number of stations with
coordinates (TransLon, TransLat).  The coordinates of my data are
(Longitude, Latitude).  Of the latter I have many more points, so the
two datasets do not share a common shape, so Ryo's method will not work
in my case.

The allowed domain is the union of subsets of the neighbourhoods of the
transect station coordinates.  So for a transect that is monotonic in
latitude, this could be done as follows:

let TransLon = Longitude[d=transect.csv]
let TransLat = Latitude[d=transect.csv]
let dx = 1
let dy = 3
set data many_observations_throughout_the_ocean.csv
let nStations = `TransLon,return=size`
def symb TransMask = (Longitude ge TransLon[i=1]-dx and Longitude le TransLon[i=1]+dx\
                  and Latitude  ge TransLat[i=1]-dy and Latitude  le TransLat[i=1]+dy)
repeat/range=2:`nStations`/name=ind ( \
    def symb TransMask = ($TransMask) or \
             (Longitude ge TransLon[i=`ind`]-dx and Longitude le TransLon[i=`ind`]+dx\
          and Latitude  ge TransLat[i=`ind`]-dy and Latitude  le TransLat[i=`ind`]+dy))
def symb myvar = if ($TransMask) then ($myvar) else 1/0 
go polymark polygon/over Latitude, Depth, "($myvar)" circle 1.0

Typically one would choose the variance of the monotonic coordinate
(dy in °) somewhat higher than the other one (dx in °), as you want to
be sure any intermediate points that are more or less on the transect
are plotted, but points far from the transect should not be plotted.

I initialise the definition with the condition for the first set of
station coordinates.  Then I add subsequent coordinates with a loop
starting from 2 upto the number of stations, by appending "OR ..." to
the condition.  I redefine the myvar symbol and plot that.  The
citation marks around ($myvar) are necessary because the symbol
contains spaces and otherwise the POLYGON subcommand only sees the
first IF instead of the whole ($TransMask) criterium as the third and
final argument (interpreting circle 1.0 as POLYMARK arguments).

There is still a problem.  With more than 12 stations it appears that I
go beyond the maximum allowed symbol length:

 **ERROR: invalid command: symbol substitution makes line too long

I can increase this by several stations by removing whitespace in my
script, but typically I have more than fifty stations that I would like
to plot.  What does it mean that the line is too long?  I can SAY a
very large string e.g. without any problems.

Marco