source: lmdz_wrf/trunk/tools/nautical.py @ 2828

## Python sript to generate nautical content
# L. Fita, CIMA. June 2019
# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
#
# pyNCplot and its component geometry_tools.py comes with ABSOLUTELY NO WARRANTY. 
# This work is licendes under a Creative Commons 
#   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
#
import os
import generic_tools as gen
import geometry_tools as geo
import numpy as np
import numpy.ma as ma
import module_ForSci as fsci

errormsg = 'ERROR -- error -- ERROR -- error'
infmsg = 'INFORMATION -- information -- INFORMATION -- information'

## Shapes/objects
# boatnames: Function to provide the names of the sections of a boat
# buoy1: Function to draw a buoy as superposition of prism and section of ball
# band_lighthouse: Function to plot a lighthouse with spiral bands
# EstuarioRioPlata: Function to plot an eschematic representation of the Estuario of Rio de la Plata
# green_buoy1: Function to draw a green mark buoy using buoy1
# isolateddanger_buoy1: Function to draw an isolated danger buoy using buoy1
# prefchannelport[A/B]_buoy1: Function to draw a preferred channel port system 
#   [A/B] buoy using buoy1
# prefchannelstarboard[A/B]_buoy1: Function to draw a preferred channel starboard 
#   system [A/B] buoy using buoy1
# red_buoy1: Function to draw a red mark buoy using buoy1
# safewater_buoy1: Function to draw a safe water mark buoy using buoy1
# special_buoy1: Function to draw an special mark buoy using buoy1
# yboat: Function to define an schematic boat from the y-plane
# z_boat: Function to define an schematic boat from the z-plane
# zsailing_boat: Function to define an schematic sailing boat from the z-plane with sails
# zisland1: Function to draw an island from z-axis as the union of a series of points by
#   circular segments
# [north/east/south/west_buoy1: Function to draw a [North/East/South/West] danger buoy using buoy1

# Definitions [Name, lat, lon]
NotablePoints = {                                                                    \
  'ArroyoRosario': ['Arroyo Rosario', np.array([-34.4331, -57.3504])],               \
  'BsAs': ['Buenos Aires', np.array([-34.6097, -58.4494])],                          \
  'BaSamborombom': ['Bah' + unichr(237) + 'a Samboromb' + unichr(243) + 'm',         \
     np.array([-36.0, -57.])],                                                       \
  'CaboPolonio': ['Cabo Polonio', np.array([-34.4083, -53.7782])],                   \
  'Colonia': ['Colonia del Sacramento', np.array([-34.4724, -57.8556])],             \
  'MartinChico': ['Martin Chico', np.array([-34.1681, -58.2118])],                   \
  'Montevideo': ['Montevideo', np.array([-34.9216, -56.1574])],                      \
  'PuntaAtalaya': ['Punta Atalaya', np.array([-35.01868, -57.5181])],                \
  'PuntaEste': ['Punta del Este', np.array([-34.9830, -54.9533])],                   \
  'PuntaIndio': ['Punta Indio', np.array([-35.4179, -57.0959])],                     \
  'PuntaMedanos': ['Punta Medanos', np.array([-36.8494, -56.6395])],                 \
  'PuntaRaza': ['Punta Raza', np.array([-36.2929, -56.7474])],                       \
  'RioSalado': ['Rio Salado', np.array([-35.7423, -57.3635])],                       \
  'Tigre': ['Tigre', np.array([-34.4486, -58.4989])],                                \
  }

# FROM: http://www.photographers1.com/Sailing/NauticalTerms&Nomenclature.html
def yboat(length=10., fcab=0.3, hcab=0.5, flength=0.7, freeboard=2., hskeg=2.,       \
  fskeg=0.2, N=200):
    """ Function to define an schematic boat from the y-plane
      length: length of the boat (without stern, default, 10)
      fcab: length of the cabin as percentage of length (default, 0.3)
      hcab: height of the cabin (default, 0.5)
      flength: floating length of the boat as percentage of length (defatult, 0.7)
      freeboard: height above the water (default, 2)
      hskeg: height of the skeg (default, 2)
      fskeg: length of the skeg as percentage of length (default, 0.2)
      N: number of points to use (default, 200)
    """
    fname = 'yboat'

    lflength = length*flength
    ilf3 = length*(1.-flength)/3
    lcab = length*fcab
    ilcab3 = length*(1.-fcab)/4.

    bow = np.array([length, freeboard])
    hbow = np.array([lflength + ilf3, 0.])
    icab = np.array([ilcab3, freeboard])
    ihcab = np.array([ilcab3, freeboard+hcab])
    ecab = np.array([ilcab3+lcab, freeboard])
    sternp = np.array([0., freeboard])
    sternlp = np.array([0., freeboard*0.8])

    print 'hbow', hbow, 'bow', bow, 'icab', icab, 'ihcab', ihcab, 'ecab', ecab, 'sternp', sternp, 'sternlp', sternlp

    boat = np.zeros((N,2), dtype=np.float)
    N1 = int(N*0.8)
    N2 = N - N1 - 1

    # stern
    N14 = N1/4
    stern = np.zeros((N14,2), dtype=np.float)
    ipt = sternlp
    ept = sternp
    dy = (ept[0] - ipt[0])/(N14-1)
    dz = (ept[1] - ipt[1])/(N14-1)
    for ip in range(N14):
        stern[ip,:] = ipt + [dy*ip, dz*ip]

    # deck
    deck = np.zeros((N14,2), dtype=np.float)
    N144 = int(N14/4.)

    ipt = sternp
    ept = icab
    dy = (ept[0] - ipt[0])/(N144-1)
    dz = (ept[1] - ipt[1])/(N144-1)
    for ip in range(N144):
        deck[ip,:] = ipt + [dy*ip, dz*ip]
    ipt = icab
    ept = ihcab
    dy = (ept[0] - ipt[0])/(N144-1)
    dz = (ept[1] - ipt[1])/(N144-1)
    for ip in range(N144):
        deck[N144:2*N144,:] = ipt + [dy*ip, dz*ip]
    deck[2*N144:3*N144,:] = geo.circ_sec(ihcab, ecab, 2*lcab, arc='short',           \
      pos='right', Nang=N144)
    N1442 = N14 - 3*N144
    ipt = ecab
    ept = bow
    dy = (ept[0] - ipt[0])/(N144-1)
    dz = (ept[1] - ipt[1])/(N144-1)
    for ip in range(N144):
        deck[3*N144:N14,:] = ipt + [dy*ip, dz*ip]

    # sternl
    sternl = np.zeros((N14,2), dtype=np.float)
    ipt = bow
    ept = hbow
    dy = (ept[0] - ipt[0])/(N14-1)
    dz = (ept[1] - ipt[1])/(N14-1)
    for ip in range(N14):
        sternl[ip,:] = ipt + [dy*ip, dz*ip]

    # keel
    N12 = N1 - 3*N14
    keel = geo.circ_sec(hbow, sternlp, length, arc='short', pos='right', Nang=N12)

    # skeg
    lskeg = length*fskeg
    ilk3 = length*(1.-fskeg)/3
    iuskeg = np.array([1.5*ilk3, 0.])
    euskeg = np.array([1.5*ilk3+lskeg, 0.])
    edskeg = np.array([1.5*ilk3+lskeg*0.8, -hskeg])
    idskeg = np.array([1.5*ilk3+lskeg*0.3, -hskeg])

    skeg = np.zeros((N2,2), dtype=np.float)
    N24=N2/4

    # upper skeg
    uskeg = np.zeros((N24,2), dtype=np.float)
    ipt = iuskeg
    ept = euskeg
    dy = (ept[0] - ipt[0])/(N24-1)
    dz = (ept[1] - ipt[1])/(N24-1)
    for ip in range(N24):
        uskeg[ip,:] = ipt + [dy*ip, dz*ip]

    # aft skeg
    askeg = np.zeros((N24,2), dtype=np.float)
    ipt = euskeg
    ept = edskeg
    dy = (ept[0] - ipt[0])/(N24-1)
    dz = (ept[1] - ipt[1])/(N24-1)
    for ip in range(N24):
        askeg[ip,:] = ipt + [dy*ip, dz*ip]

    # down skeg
    dskeg = np.zeros((N24,2), dtype=np.float)
    ipt = edskeg
    ept = idskeg
    dy = (ept[0] - ipt[0])/(N24-1)
    dz = (ept[1] - ipt[1])/(N24-1)
    for ip in range(N24):
        dskeg[ip,:] = ipt + [dy*ip, dz*ip]

    # stern skeg
    N22 = N2 - 3*N24
    sskeg = np.zeros((N22,2), dtype=np.float)
    ipt = idskeg
    ept = iuskeg
    dy = (ept[0] - ipt[0])/(N22-1)
    dz = (ept[1] - ipt[1])/(N22-1)
    for ip in range(N22):
        sskeg[ip,:] = ipt + [dy*ip, dz*ip]

    boat[0:N14,:] = stern
    boat[N14:2*N14,:] = deck
    boat[2*N14:3*N14,:] = sternl
    boat[3*N14:4*N12,:] = keel
    boat[N1,:] = np.array([gen.fillValueF, gen.fillValueF])
    boat[N1+1:N1+1+N24,:] = uskeg
    boat[N1+1+N24:N1+1+2*N24,:] = askeg
    boat[N1+1+2*N24:N1+1+3*N24,:] = dskeg
    boat[N1+1+3*N24:N,:] = sskeg

    # correct order of sections
    boatsecs = ['stern', 'deck', 'sternl', 'keel', 'uskeg', 'askeg', 'dskeg', 'sskeg']

    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
    dicboat = {'stern': [stern, '-', '#8A5900', 2.],                         \
      'deck': [deck, '-', '#8A5900', 2.],                                  \
      'sternl': [sternl, '-', '#8A5900', 2.],                                          \
      'keel': [keel, '-', '#8A5900', 2.],                        \
      'uskeg': [uskeg, '-', '#000000', 1.5], 'askeg': [askeg, '-.', '#000000', 1.5], \
      'dskeg': [dskeg, '-', '#000000', 1.5], 'sskeg': [sskeg, '-.', '#000000', 1.5]}

    boat = ma.masked_equal(boat, gen.fillValueF)
      
    return boat, boatsecs, dicboat

def zboat(length=10., beam=1., lbeam=0.4, sternbp=0.5):
    """ Function to define an schematic boat from the z-plane
      length: length of the boat (without stern, default 10)
      beam: beam of the boat (default 1)
      lbeam: length at beam (as percentage of length, default 0.4)
      sternbp: beam at stern (as percentage of beam, default 0.5)
    """
    fname = 'zboat'

    bow = np.array([length, 0.])
    maxportside = np.array([length*lbeam, -beam])
    maxstarboardside = np.array([length*lbeam, beam])
    portside = np.array([0., -beam*sternbp])
    starboardside = np.array([0., beam*sternbp])

    # forward section
    fportside = geo.circ_sec(maxportside, bow, length*2)
    fstarboardside = geo.circ_sec(bow, maxstarboardside, length*2)
    # aft section
    aportside = geo.circ_sec(portside, maxportside, length*2)
    astarboardside = geo.circ_sec(maxstarboardside, starboardside, length*2)
    # stern
    stern = geo.circ_sec(starboardside, portside, length*2)

    dpts = stern.shape[0]
    boat = np.zeros((dpts*5,2), dtype=np.float)

    boat[0:dpts,:] = aportside
    boat[dpts:2*dpts,:] = fportside
    boat[2*dpts:3*dpts,:] = fstarboardside
    boat[3*dpts:4*dpts,:] = astarboardside
    boat[4*dpts:5*dpts,:] = stern

    fname = 'boat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) + '_lb' +  \
      str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) + '.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
          ' %\n')
        for ip in range(dpts*5):
            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"


    # Center line extending [fcl] percentage from length on aft and stern
    fcl = 0.15
    centerline = np.zeros((dpts,2), dtype=np.float)
    dl = length*(1.+fcl*2.)/(dpts-1)
    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)

    # correct order of sections
    boatsecs = ['aportside', 'fportside', 'fstarboardside', 'astarboardside',        \
      'stern', 'centerline']

    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
    dicboat = {'fportside': [fportside, '-', '#8A5900', 2.],                         \
      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
      'stern': [stern, '-', '#8A5900', 2.],                                          \
      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
      'centerline': [centerline, '-.', '#AA6464', 1.5]}
    
    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +'.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' +str(sternbp)+'\n')
        for ip in range(dpts*5):
            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"
  
    return boat, boatsecs, dicboat

def zsailing_boat(length=10., beam=1., lbeam=0.4, sternbp=0.5, lmast=0.6, wmast=0.1, \
  hsd=5., msd=5., lheads=0.38, lmains=0.55):
    """ Function to define an schematic sailing boat from the z-plane with sails
      length: length of the boat (without stern, default 10)
      beam: beam of the boat (default 1)
      lbeam: length at beam (as percentage of length, default 0.4)
      sternbp: beam at stern (as percentage of beam, default 0.5)
      lmast: position of the mast (as percentage of length, default 0.6)
      wmast: width of the mast (default 0.1)
      hsd: head sail direction respect to center line (default 5., -999.99 for upwind)
      msd: main sail direction respect to center line (default 5., -999.99 for upwind)
      lheads: length of head sail (as percentage of legnth, defaul 0.38)
      lmains: length of main sail (as percentage of legnth, defaul 0.55)
    """
    fname = 'zsailing_boat'

    bow = np.array([length, 0.])
    maxportside = np.array([length*lbeam, -beam])
    maxstarboardside = np.array([length*lbeam, beam])
    portside = np.array([0., -beam*sternbp])
    starboardside = np.array([0., beam*sternbp])

    aportside = geo.circ_sec(portside, maxportside, length*2)
    fportside = geo.circ_sec(maxportside, bow, length*2)
    fstarboardside = geo.circ_sec(bow, maxstarboardside, length*2)
    astarboardside = geo.circ_sec(maxstarboardside, starboardside, length*2)
    stern = geo.circ_sec(starboardside, portside, length*2)
    dpts = fportside.shape[0]

    # correct order of sections
    sailingboatsecs = ['aportside', 'fportside', 'fstarboardside', 'astarboardside', \
      'stern', 'mast', 'hsail', 'msail', 'centerline']

    # forward section

    # aft section
    # stern
    # mast
    mast = geo.p_circle(wmast,N=dpts)
    mast = mast + [length*lmast, 0.]
    # head sails
    lsail = lheads*length
    if hsd != -999.99:
        sailsa = np.pi/2. - np.pi*hsd/180.
        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
        endsail[0] = length - endsail[0]
        if bow[1] > endsail[1]:
            hsail = geo.circ_sec(endsail, bow, lsail*2.15)
        else:
            hsail = geo.circ_sec(bow, endsail, lsail*2.15)
    else:
        hsail0, sailsec, saildic = geo.p_sinusiode(length=lsail, amp=0.2, lamb=0.75, N=dpts)
        hsail = np.zeros((dpts,2), dtype=np.float)
        hsail[:,0] = hsail0[:,1]
        hsail[:,1] = hsail0[:,0]
        hsail = bow - hsail

    # main sails
    lsail = lmains*length
    if msd != -999.99:
        sailsa = np.pi/2. - np.pi*msd/180.
        begsail = np.array([length*lmast, 0.])
        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
        endsail[0] = length*lmast - endsail[0]
        if endsail[1] > begsail[1]:
            msail = geo.circ_sec(begsail, endsail, lsail*2.15)
        else:
            msail = geo.circ_sec(endsail, begsail, lsail*2.15)
    else:
        msail0, sailsec, saildic = geo.p_sinusiode(length=lsail, amp=0.25, lamb=1., N=dpts)
        msail = np.zeros((dpts,2), dtype=np.float)
        msail[:,0] = msail0[:,1]
        msail[:,1] = msail0[:,0]
        msail = [length*lmast,0] - msail

    sailingboat = np.zeros((dpts*8+4,2), dtype=np.float)

    sailingboat[0:dpts,:] = aportside
    sailingboat[dpts:2*dpts,:] = fportside
    sailingboat[2*dpts:3*dpts,:] = fstarboardside
    sailingboat[3*dpts:4*dpts,:] = astarboardside
    sailingboat[4*dpts:5*dpts,:] = stern
    sailingboat[5*dpts,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[5*dpts+1:6*dpts+1,:] = mast
    sailingboat[6*dpts+1,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[6*dpts+2:7*dpts+2,:] = hsail
    sailingboat[7*dpts+2,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[7*dpts+3:8*dpts+3,:] = msail
    sailingboat[8*dpts+3,:] = [gen.fillValueF, gen.fillValueF]

    sailingboat = ma.masked_equal(sailingboat, gen.fillValueF)

    # Center line extending [fcl] percentage from length on aft and stern
    fcl = 0.15
    centerline = np.zeros((dpts,2), dtype=np.float)
    dl = length*(1.+fcl*2.)/(dpts-1)
    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)

    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
    dicsailingboat = {'fportside': [fportside, '-', '#8A5900', 2.],                  \
      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
      'stern': [stern, '-', '#8A5900', 2.],                                          \
      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
      'mast': [mast, '-', '#8A5900', 2.], 'hsail': [hsail, '-', '#AAAAAA', 1.],      \
      'msail': [msail, '-', '#AAAAAA', 1.],                                          \
      'centerline': [centerline, '-.', '#AA6464', 1.5]}
    
    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +                \
      '_lm' + str(int(lmast*100.)) + '_wm' + str(int(wmast)) +                       \
      '_hsd' + str(int(hsd)) + '_hs' + str(int(lheads*100.)) +                       \
      '_ms' + str(int(lheads*100.)) + '_msd' + str(int(msd)) +'.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
          ' % mast position: '+ str(lmast) + ' % mast width: ' + str(wmast) + ' ' +  \
          ' head sail direction:' + str(hsd) + ' head sail length: ' + str(lheads) + \
          ' %' + ' main sail length' + str(lmains) + ' main sail direction:' +       \
          str(msd) +'\n')
        for ip in range(dpts*5):
            of.write(str(sailingboat[ip,0]) + ' ' + str(sailingboat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"
 
    return sailingboat, sailingboatsecs, dicsailingboat

def zisland1(mainpts= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],\
  [2.8, -0.1], [0.1,-0.6]], dtype=np.float), radfrac=3., N=200):
    """ Function to draw an island from z-axis as the union of a series of points by
        circular segments
      mainpts: main points of the island (clockwise ordered, to be joined by 
        circular segments of radii as the radfrac factor of the distance between 
        consecutive points)
          * default= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],
            [2.8, -0.1], [0.1,-0.6]], dtype=np.float)
      radfrac: multiplicative factor of the distance between consecutive points to 
        draw the circular segment (3., default)
      N: number of points (200, default)
    """
    fname = 'zisland1'

    island1 = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # Coastline
    island1 = geo.join_circ_sec_rand(mainpts, arc='short', pos='left')

    islandsecs = ['coastline']
    islanddic = {'coastline': [island1, '-', '#161616', 2.]}

    island1 = ma.masked_equal(island1, gen.fillValueF)

    return island1, islandsecs, islanddic

def buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=300):
    """ Function to draw a buoy as superposition of prism and section of ball
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'buoy1'

    buoy = np.zeros((N,2), dtype=np.float)

    N3 = int(N/3/5)
    NNp = 0
    iip = 0
    # left lateral
    ix = -width/2.
    Np = N3
    iy = 0.
    dx = 0.
    dy = height/(Np)
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # left upper
    ix = -width/2.
    iy = height
    dx = (width/2.-bradii*bfrac)/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # ball
    p1 = np.array([height, -bradii*bfrac])
    p2 = np.array([height, bradii*bfrac])
    Np = int(2*N/3)
    buoy[iip:iip+Np,:] = geo.circ_sec(p1, p2, 2.*bradii, 'long', 'left', Np)
    NNp = NNp + Np
    iip = NNp

    # right upper
    ix = bradii*bfrac
    iy = height
    Np = N3
    dx = (width/2.-bradii*bfrac)/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # right lateral
    ix = width/2.
    iy = height
    dx = 0.
    dy = -height/(Np)
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # Base
    ix = width/2.
    iy = 0.
    Np = N - int(2*N/3) - 4*N3 - 1
    dx = -width/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    buoy[N-1,:] = buoy[0,:]

    buoysecs = ['base']
    buoydic = {'base': [buoy, '-', 'k', 1.5]}

    return buoy, buoysecs, buoydic

def band_lighthouse(height=10., width=2., hlight=3., bands=3, N=300):
    """ Function to plot a lighthouse with spiral bands
      height: height of the tower (10., default)
      width: width of the tower (2., default)
      hlight: height of the light (3., default)
      bands: number of spiral bands (3, default)
      N: number of points (300, default)
    """
    fname = 'band_lighthouse'

    lighthouse = np.ones((N,2), dtype=np.float)*gen.fillValueF
    lighthousesecs = []
    lighthousedic = {}

    # base Tower
    Nsec = int(0.30*N/7)
    p1=np.array([0., width/2.])
    p2=np.array([0., -width/2.])
    iip = 0
    lighthouse[0:Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec

    # left side
    ix=-width/2.
    iy=0.
    dx = 0.
    dy = height/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec

    # Top Tower
    p1=np.array([height, width/2.])
    p2=np.array([height, -width/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec

    # right side
    ix=width/2.
    iy=height
    dx = 0.
    dy = -height/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec + 1

    Ntower = iip-1
    lighthousesecs.append('tower')
    lighthousedic['tower'] = [lighthouse[0:iip-1], '-', 'k', 1.5]

    # Left light
    p1 = np.array([height, -width*0.8/2.])
    p2 = np.array([height+hlight, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*hlight, Nang=Nsec)
    iip = iip + Nsec
    
    # Top Light
    p1=np.array([height+hlight, width*0.8/2.])
    p2=np.array([height+hlight, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec + 1

    # Right light
    p1 = np.array([height+hlight, width*0.8/2.])
    p2 = np.array([height, width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*hlight, Nang=Nsec)
    iip = iip + Nsec

    # Base Light
    p1=np.array([height, width*0.8/2.])
    p2=np.array([height, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec + 1
    lighthousesecs.append('light')
    lighthousedic['light'] = [lighthouse[Ntower+1:iip-1], '-', '#EEEE00', 1.5]

    # Spiral bands
    hb = height/(2.*bands)
    Nsec2 = (N - Nsec*8 - 3)/bands
    for ib in range(bands-1):
        iband = iip
        Nsec = Nsec2/4
        bandS = 'band' + str(ib).zfill(2)
        # hband
        ix = -width/2.
        iy = hb*ib*2
        dx = 0.
        dy = hb/(Nsec-1)
        for ip in range(Nsec):
            lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
        iip = iip + Nsec
        # uband
        p1 = np.array([hb*(ib*2+1), -width/2.])
        p2 = np.array([hb*(ib*2+2), width/2.])
        lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='right', Nang=Nsec)
        iip = iip + Nsec
        # dband
        ix = width/2.
        iy = hb*(ib*2+2)
        dx = 0.
        dy = -hb/(Nsec-1)
        for ip in range(Nsec):
            lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
        iip = iip + Nsec
        # dband
        p1 = np.array([hb*(ib*2+1), width/2.])
        p2 = np.array([hb*ib*2, -width/2.])
        lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
        iip = iip + Nsec + 1
        lighthousesecs.append(bandS)
        lighthousedic[bandS] = [lighthouse[iband:iip-1], '-', '#6408AA', 2.]

    ib = bands-1
    Nsec3 = (N - iip - 1)
    Nsec = int(Nsec3/4)
    bandS = 'band' + str(ib).zfill(2)
    # hband
    iband = iip
    ix = -width/2.
    iy = hb*ib*2
    dx = 0.
    dy = hb/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec
    # uband
    p1 = np.array([hb*(ib*2+1), -width/2.])
    p2 = np.array([hb*(ib*2+2), width/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='right', Nang=Nsec)
    iip = iip + Nsec
    # dband
    ix = width/2.
    iy = hb*(2+ib*2)
    dx = 0.
    dy = -hb/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec
    # dband
    Nsec = N - iip
    p1 = np.array([hb*(1+ib*2), width/2.])
    p2 = np.array([hb*ib*2, -width/2.])
    lighthouse[iip:iip+Nsec,:] = geo.circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)        
    lighthousesecs.append(bandS)
    lighthousedic[bandS] = [lighthouse[iband:iip-1], '-', '#6408AA', 2.]

    lighthouse = ma.masked_equal(lighthouse, gen.fillValueF)

    return lighthouse, lighthousesecs, lighthousedic

def north_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a North danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'north_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    triu = geo.p_angle_triangle(N=N32)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.1*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    trid = geo.p_angle_triangle(N=N323)
    trib = trid*lsign + [0.,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+1.1*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = geo.cut_ypolygon(buoy1v, yval=Height/2., keep='below')
    Ncut, halfup = geo.cut_ypolygon(buoy1v, yval=Height/2., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign1', 'sign2', 'half1', 'half2']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'half1': [halfup, '-', 'k', 1.],    \
      'half2': [halfdown, '-', '#FFFF00', 1.]}

    return buoy, buoysecs, buoydic

def east_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a East danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'east_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    triu = geo.p_angle_triangle(N=N32)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.1*lsign,0.]

    # down
    N323 = N - N32 - N2 - 2

    trid = geo.p_angle_triangle(N=N323)
    trid = geo.mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+0.9*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, halfbtw = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, halfup = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign1', 'sign2', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5],                                     \
      'third1': [halfup, '-', 'k', 1.], 'third2': [halfbtw, '-', '#FFFF00', 1.],     \
      'third3': [halfdown, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def south_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a South danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'south_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    trid = geo.p_angle_triangle(N=N32)
    trid = geo.mirror_polygon(trid, 'x')
    trib = trid*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.9*lsign,0.]

    # down
    N323 = N - N32 - N2 - 2
    trid = geo.p_angle_triangle(N=N323)
    trid = geo.mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+0.9*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = geo.cut_ypolygon(buoy1v, yval=Height/2., keep='below')
    Ncut, halfup = geo.cut_ypolygon(buoy1v, yval=Height/2., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign1', 'sign2', 'half1', 'half2']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'half1': [halfup, '-', '#FFFF00', 1.], \
      'half2': [halfdown, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def west_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a West danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'east_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns

    # down
    N32 = int(N3/2) 
    trid = geo.p_angle_triangle(N=N32)
    trid = geo.mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+1.9*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    triu = geo.p_angle_triangle(N=N323)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+N323+2:N,:] = trib + [bottsigns+1.*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, halfbtw1 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, halfup = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign1', 'sign2', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'third1': [halfdown, '-', '#FFFF00', 1.], 'third2': [halfbtw1, '-', 'k', 1.],  \
      'third3': [halfup, '-', '#FFFF00', 1.],                                        \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5]}

    return buoy, buoysecs, buoydic

def safewater_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a safe water mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'safewater_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns
    
    Height = np.max(buoy1v[:,0])
    sign = geo.p_circle(lsign, N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    ix = -width/2.
    Ncut, quarter1 = geo.cut_xpolygon(buoy1v, xval=ix+width/4., keep='left')
    Ncut, quarter2 = geo.cut_between_xpolygon(buoy1v, xval1=ix+width/4., xval2=ix+width/2.)
    Ncut, quarter3 = geo.cut_between_xpolygon(buoy1v, xval1=ix+width/2., xval2=ix+3.*width/4.)
    Ncut, quarter4 = geo.cut_xpolygon(buoy1v, xval=ix+3.*width/4., keep='right')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'quarter1', 'quarter2', 'quarter3', 'quarter4']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'quarter1': [quarter1,'-','r',1.], \
      'quarter2': [quarter2,'-','#FFFFFF',1.], 'quarter3': [quarter3,'-','r',1.],    \
      'quarter4': [quarter4,'-','#FFFFFF',1.]}

    return buoy, buoysecs, buoydic

def red_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a red mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'red_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = geo.p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5]}

    return buoy, buoysecs, buoydic

def green_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a green mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'green_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = geo.p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelportA_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, \
  N=300):
    """ Function to draw a preferred channel port system A buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelportA_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = geo.p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    Ncut, third1 = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5], 'third1': [third1,'-','g',1.5],    \
      'third2': [third2,'-','r',1.5], 'third3': [third3,'-','g',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelportB_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, \
  N=300):
    """ Function to draw a preferred channel port system B buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelportB_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = geo.p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    Ncut, third1 = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'third1': [third1,'-','r',1.5],    \
      'third2': [third2,'-','g',1.5], 'third3': [third3,'-','r',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelstarboardA_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8,        \
  hsigns=0.3, N=300):
    """ Function to draw a preferred channel starboard system A buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelstarboardA_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = geo.p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    # painting it
    Ncut, third1 = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'third1': [third1,'-','r',1.5],    \
      'third2': [third2,'-','g',1.5], 'third3': [third3,'-','r',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelstarboardB_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8,        \
  hsigns=0.3, N=300):
    """ Function to draw a preferred channel starboard system B buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelstarboardB_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = geo.p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    # painting it
    Ncut, third1 = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5], 'third1': [third1,'-','g',1.5],    \
      'third2': [third2,'-','r',1.5], 'third3': [third3,'-','g',1.5]}

    return buoy, buoysecs, buoydic

def isolateddanger_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5,   \
  N=300):
    """ Function to draw an isolated danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'isolateddanger_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    circle = geo.p_circle(lsign/2., N=N32)
    trib = circle + [0.,0.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+3.2*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    trid = geo.p_circle(lsign/2., N=N32)
    trib = circle + [0.,0.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+2.*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, third1 = geo.cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = geo.cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = geo.cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign1', 'sign2', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'third1': [third1, '-', 'k', 1.],   \
      'third2': [third2, '-', 'r', 1.], 'third3': [third3, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def special_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5, N=300):
    """ Function to draw an special mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'special_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    Height = np.max(buoy1v[:,0])

    # sign
    N3 = N - N2 - 1
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    cross, crosssecs, crossdic = geo.p_cross_width(lsign, width=0.3*lsign, Narms=2, N=N3)
    cross = geo.rotate_polygon_2D(cross, 40.05)
    buoy[N2+1:N,:] = cross + [Height+1.1*lsign,0.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','#FFFF00',1.5],                             \
      'sign': [buoy[N2+1:N,:],'-','#FFFF00',1.5]}

    return buoy, buoysecs, buoydic

def emergency_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5, N=300):
    """ Function to draw an eergency mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'emergency_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    Height = np.max(buoy1v[:,0])

    # sign
    N3 = N - N2 - 1
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    cross, crosssecs, crossdic = geo.p_cross_width(lsign, width=0.3*lsign, Narms=2, N=N3)
    buoy[N2+1:N,:] = cross + [Height+1.1*lsign,0.]

    # painting it
    ix = -width/2.
    Ncut, fifth1 = geo.cut_xpolygon(buoy1v, xval=ix+width/5., keep='left')
    Ncut, fifth2 = geo.cut_between_xpolygon(buoy1v,xval1=ix+width/5.,xval2=ix+width*2./5.)
    Ncut, fifth3 = geo.cut_between_xpolygon(buoy1v,xval1=ix+width*2./5.,xval2=ix+width*3./5.)
    Ncut, fifth4 = geo.cut_between_xpolygon(buoy1v,xval1=ix+width*3./5.,xval2=ix+width*4./5.)
    Ncut, fifth5 = geo.cut_xpolygon(buoy1v, xval=ix+width*4./5., keep='right')

    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign', 'fifth1', 'fifth2', 'fifth3', 'fifth4', 'fifth5']
    buoydic = {'buoy': [buoy[0:N2,:],'-','#FFFF00',1.5],                             \
      'sign': [buoy[N2+1:N,:],'-','#FFFF00',1.5],'fifth1':[fifth1,'-','#FFFF00',1.5],\
      'fifth2': [fifth2,'-','#0000FF',1.5],'fifth3': [fifth3,'-','#FFFF00',1.5],     \
      'fifth4': [fifth4,'-','#0000FF',1.5],'fifth5': [fifth5,'-','#FFFF00',1.5]}

    return buoy, buoysecs, buoydic

def EstuarioRioPlata(N=300):
    """ Function to plot an eschematic representation of the Estuario of Rio de la Plata
      N: total number of vertices to use
    """
    fname = 'EstuarioRioPlata'

    secs0 = ['PuntaMedanos', 'PuntaRaza', 'RioSalado', 'PuntaIndio', 'PuntaAtalaya', \
      'Tigre', 'MartinChico', 'Colonia', 'ArroyoRosario', 'Montevideo', 'PuntaEste', \
      'CaboPolonio']
    secs = []
    dic = {}
    rads = [5., 1.0, 5., 5., 5., 5., 5., 5., 5., 5., 5.]
    lengths = ['short', 'short', 'short', 'short', 'short', 'short', 'short',        \
      'short', 'short', 'short', 'short']
    sides = ['right', 'left', 'left', 'right', 'left', 'left', 'left', 'left',       \
      'right', 'left', 'right']
    Nsecs = len(secs0)
    Nn = N/Nsecs
    estuario = np.zeros((N,2), dtype=np.float)

    iip = 0
    # Atlantic_PuntaRaza
    prevn = 'PuntaMedanos'
    pv = NotablePoints[prevn]
    ip = pv[1]
    for isec in range(1,Nsecs-1):
        iisec = isec - 1
        aname = secs0[isec]
        pv = NotablePoints[aname]
        ep = pv[1]
        dps = geo.dist_points(ip,ep)
        estuario[iip:iip+Nn,:] = geo.circ_sec(ip,ep, dps*rads[iisec], lengths[iisec],\
          sides[iisec], Nn)
        secs.append(prevn+'_'+aname)
        dic[prevn+'_'+aname] = [estuario[iip:iip+Nn,:], ['-', 'k', 1.]]
        ip = ep + 0.
        prevn = aname + ''
        iip = iip + Nn

    Nn2 =  N - (Nsecs-2)*Nn
    isec = Nsecs-1
    iisec = isec - 1
    aname = secs0[isec]
    pv = NotablePoints[aname]
    ep = pv[1]
    dps = geo.dist_points(ip,ep)
    isec = Nsecs - 1
    estuario[iip:N,:] = geo.circ_sec(ip, ep, dps*rads[iisec], lengths[iisec],        \
      sides[iisec], Nn2)
    secs.append(prevn+'_'+aname)
    dic[prevn+'_'+aname] = [estuario[iip:N,:], ['-', 'k', 1.]]

    return estuario, secs, dic

def boatnames(xn,xx,yn,yx,zn,zx,zlf):
    """ Function to provide the names of the sections of a boat
      xn: minimum length on x-axis (across beam)
      xx: maximum length on x-axis (across beam)
      yn: minimum length on y-axis (length)
      yx: maximum length on y-axis (length)
      zn: minimum length on z-axis (draught)
      zx: maximum length on z-axis (draught)
      zlf: water line
    """
    fname = 'boatnames'

    dx = xx - xn
    dy = yx - yn
    dz = zx - zn

    x0 = xn + dx/2.
    y0 = yn + dy/2.
    z0 = zn + dz/2.

    # Values
    boatvs = {
      'xn': xn, 'xx': xx, 'yn': yn, 'yx': yx, 'zn': zn, 'zx': zx,                    \
      'dx': dx, 'dy': dy, 'dz': dz, 'zlf': zlf,                                      \
      }

    # Names
    boatns = {
      'bow': ['bow', 'proa', np.array([x0,yx,zx])],                                  \
      'stern': ['stern', 'popa', np.array([x0,yn,zx])],                              \
      'starboard': ['starboard', 'estribor', np.array([xx,y0,zx])],                  \
      'port': ['port', 'babor', np.array([xn,y0,zx])],                               \
      'waterline': ['waterline', 'l'+unichr(237)+'nea de flotaci'+ unichr(243)+'n',  \
        np.array([xn,y0,zlf])],                                                      \
      'keel': ['keel', 'quillote', np.array([xn,y0,zn])],                            \
      'centerline': ['center line', 'l'+unichr(237)+'nea de cruj'+unichr(237)+       \
        'a (plano)', np.array([x0,y0,zn])],                                         \
      'bowside': ['bow', 'amura',  np.array([xx,yx*0.83,zx])],                       \
      'beamside': ['beam', 'trav' + unichr(233)+ 's',  np.array([xx,yx*0.5,zx])],    \
      'quarter': ['quarter', 'aleta',  np.array([xx,yx*0.15,zx])],                   \

      }

    # Dimensions
    boatls = {
      'length': ['length', 'eslora', np.array([[x0,yn,zx], [x0,yx,zx]])],            \
      'beam': ['beam', 'manga', np.array([[xn,y0,zx], [xx,y0,zx]])],                 \
      'freeboard': ['freeboard (air \ndraught)', 'francobordo (obra \nviva)\n carena', \
        np.array([[xn,yn,zlf], [xn,yn,zx]])],                                        \
      'draught': ['draught', 'calado (obra \nmuerta)',                               \
        np.array([[xn,yx,zlf],[xn,yx,zn]])],                                         \
      'bowside': ['bow', 'amura',                                                    \
        np.array([[xx,yx*0.6,zx],[xn,yx*0.6,zx]])],                                  \
      'beamside': ['beam', 'trav'+unichr(233)+'s',                                   \
        np.array([[xx,yx*0.3,zx], [xn,yx*0.3,zx]])],                                 \
      'quarter': ['quarter', 'aleta', np.array([[xx,0.,zx], [xn,0.,zx]])],           \
      }

    return boatvs, boatns, boatls