source: issm/oecreview/Archive/24684-25833/ISSM-25064-25065.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 78.5 KB

  • ../trunk-jpl/src/m/shp/shpwrite.py

     
     1try:
     2    import shapefile
     3except ImportError:
     4    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
     5
     6def shpwrite(shp, filename): #{{{
     7    '''
     8    SHPREAD - write a shape file from a contour structure
     9
     10    The current implementation of shpwrite depends on PyShp.
     11
     12    Usage:
     13        shpwrite(shp, filename)
     14
     15    Example:
     16        shpwrite(shp, 'domainoutline.shp')
     17
     18    See also SHPREAD
     19
     20    Sources:
     21    - https://github.com/GeospatialPython/pyshp
     22
     23    TODO:
     24    - Should we check if there is only one element (see how MATLAB's shaperead
     25    and shapewrite handle single shapes versus multiple shapes)?
     26    '''
     27
     28    # Sanity checks
     29    for shape in shp:
     30        print(shape)
     31
     32    if shp[0].Geometry == 'Point':
     33        shapeType = 1
     34    elif shp[0].Geometry == 'Line':
     35        shapeType = 3
     36    elif shp[0].Geometry == 'Polygon':
     37        shapeType = 5
     38    else:
     39        raise Exception('shpwrite error: geometry \'{}\' is not currently supported'.format(shp[0].Geometry))
     40
     41    sf = shapefile.Writer(filename, shapeType=shapeType)
     42
     43    for i in range(len(shp)):
     44        sf.field('name', 'C') # TODO: Allow shape ids/names to be passed through
     45        if shapeType == 1: # POINT
     46            sf.point(shp[i].x, shp[i].y)
     47        elif shapeType == 3: # POLYLINE
     48            points = []
     49            for j in range(len(shp[i].x)):
     50                points.append([shp[i].x[j], shp[i].y[j]])
     51            sf.line(line)
     52        elif shapeType == 5: # POLYGON
     53            points = []
     54            for j in range(len(shp[i].x)):
     55                points.append([shp[i].x[j], shp[i].y[j]])
     56            sf.poly(points)
     57        sf.null()
     58        sf.record(str(i))
     59    sf.close()
     60#}}}

  • ../trunk-jpl/src/m/shp/shpwrite.m

     
    99%
    1010%   See also SHPREAD
    1111
    12 
    13 %initialize number of profile
    14 count=0;
    15 
    1612contours=struct([]);
    1713for i=1:length(shp),
    1814        if strcmpi(shp(i).Geometry,'Point'),
    … …  
    2117                contours(i).Geometry='Polygon';
    2218        elseif strcmpi(shp(i).Geometry,'Line'),
    2319                contours(i).Geometry='Line';
     20        else,
     21                error(['shpwrite error: geometry ' shp(i).Geometry ' not currently supported']);
    2422        end
    2523        contours(i).id=i;
    2624        contours(i).X=shp(i).x;
    2725        contours(i).Y=shp(i).y;
    2826end
    29        
     27
    3028shapewrite(contours,filename);
    3129end

  • ../trunk-jpl/src/m/mesh/meshintersect3d.py

     
    66
    77from pairoptions import pairoptions
    88
     9
    910def meshintersect3d(x, y, z, xs, ys, zs, *args): #{{{
    1011    '''
    11     MESHINTERSECT - returns indices (into x, y, and z) of common values between (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
     12    MESHINTERSECT - returns indices (into x, y, and z) of common values between
     13    (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
    1214    '''
    1315
    1416    # Process options

  • ../trunk-jpl/src/m/coordsystems/laea.py

     
     1def laea(lat, long): #{{{
     2    '''
     3    LAEA - Lambert Azimuthal Equal Area projection at lat, long projection
     4    center.
     5
     6        Usage:
     7            string = laea(45, -90)
     8
     9        Example:
     10            string = laea(45, -90)
     11            return string = '+proj=laea +lat_0=45 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'
     12    '''
     13
     14    return '+proj=laea +lat_0={} +lon_0={} +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'.format(lat, long)
     15#}}}

  • ../trunk-jpl/src/m/coordsystems/gdaltransform.m

     
    1919%               Bamber's Greeland projection
    2020%                       +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
    2121%
    22 %               To get proj.4 string from EPSG, use gdalsrsinfo. Example:
     22%       To get proj.4 string from EPSG, use gdalsrsinfo. Example:
    2323%
    24 %                       gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
     24%               gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
    2525
    2626        %give ourselves unique file names
    2727        filename_in  = tempname();

  • ../trunk-jpl/src/m/geometry/AboveGround.py

     
     1import numpy as np
     2
     3def function AboveGround(lat, long, r, height): #{{{
     4    r = r + height
     5    x = r * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(long))
     6    y = r * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(long))
     7    z = r * np.sin(np.deg2rad(lat))
     8#}}}

  • ../trunk-jpl/src/m/classes/plotoptions.py

     
    11from collections import OrderedDict
     2
    23import pairoptions
    34
    45
    … …  
    67    '''
    78    PLOTOPTIONS class definition
    89
    9         Usage:
    10             plotoptions = plotoptions(*arg)
     10    Usage:
     11        plotoptions = plotoptions(*arg)
    1112    '''
    1213
    1314    def __init__(self, *arg):  # {{{
    … …  
    3536    def buildlist(self, *arg):  #{{{
    3637        #check length of input
    3738        if len(arg) % 2:
    38             raise TypeError('Invalid parameter / value pair arguments')
     39            raise TypeError('Invalid parameter/value pair arguments')
    3940
    4041        #go through args and build list (like pairoptions)
    4142        rawoptions = pairoptions.pairoptions(*arg)
    … …  
    102103                        if len(nums) != 2:
    103104                            continue
    104105                        if False in [x.isdigit() for x in nums]:
    105                             raise ValueError('error: in option i - j both i and j must be integers')
     106                            raise ValueError('error: in option i-j both i and j must be integers')
    106107                        for j in range(int(nums[0]) - 1, int(nums[1])):
    107108                            self.list[j].addfield(field, rawlist[i][1])
    108109

  • ../trunk-jpl/src/m/classes/basin.m

     
    1818        end
    1919        methods
    2020                function self = basin(varargin) % {{{
    21                         switch nargin
    22                                 case 0
    23                                         self=setdefaultparameters(self);
    24                                 otherwise
     21                        self=setdefaultparameters(self);
    2522
    26                                         self=setdefaultparameters(self);
    27                                         options=pairoptions(varargin{:});
    28                        
    29                                         %recover field values:
    30                                         self.boundaries=getfieldvalue(options,'boundaries',{});
    31                                         self.name=getfieldvalue(options,'name','');
    32                                         self.continent=getfieldvalue(options,'continent','');
     23                        if nargin>0,
     24                                options=pairoptions(varargin{:});
     25               
     26                                %recover field values:
     27                                self.boundaries=getfieldvalue(options,'boundaries',{});
     28                                self.name=getfieldvalue(options,'name','');
     29                                self.continent=getfieldvalue(options,'continent','');
    3330
    34                                         %figure out projection string:
     31                                %figure out projection string:
     32                                if exist(options,'epsg'),
     33                                        if exist(options,'proj'),
     34                                                error(['Error in constructor for basin ' self.name '. Cannot supply epsg and proj at the same time!']);
     35                                        end
     36                                        epsg=getfieldvalue(options,'epsg');
     37                                        proj=epsg2proj(epsg); %convert to PROJ.4 format
     38                                elseif exist(options,'proj'),
    3539                                        if exist(options,'epsg'),
    36                                                 if exist(options,'proj'),
    37                                                         error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
    38                                                 end
    39                                                 epsg=getfieldvalue(options,'epsg');
    40                                                 %convert into proj4:
    41                                                 proj=epsg2proj(epsg);
    42                                         elseif exist(options,'proj'),
    43                                                 if exist(options,'epsg'),
    44                                                         error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
    45                                                 end
    46                                                 proj=getfieldvalue(options,'proj');
    47                                         else
    48                                                 proj= '+proj=longlat +datum=WGS84 +no_defs';
     40                                                error(['Error in constructor for basin ' self.name '. Cannot supply proj and epsg at the same time!']);
    4941                                        end
    50                                         self.proj=proj;
     42                                        proj=getfieldvalue(options,'proj');
     43                                else
     44                                        proj='+proj=longlat +datum=WGS84 +no_defs';
     45                                end
     46
     47                                self.proj=proj;
    5148                        end
    5249                end % }}}
    5350                function self = setdefaultparameters(self) % {{{
    5451                        self.name='';
    5552                        self.continent='';
    56                         self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
     53                        self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326
    5754                        self.boundaries={};
    5855
    5956                end % }}}
    … …  
    9289                        options=pairoptions(varargin{:});
    9390                        extension=getfieldvalue(options,'extension',1);
    9491
    95                         [path,nme,ext]=fileparts(self.name);
     92                        [path,name,ext]=fileparts(self.name);
    9693                        if extension,
    97                                 output=[nme ext];
     94                                output=[name ext];
    9895                        else
    99                                 output=nme;
     96                                output=name;
    10097                        end
    10198                end % }}}
    10299                function plot(self,varargin) % {{{
    … …  
    173170                        inshapefile=getfieldvalue(options,'shapefile');
    174171                        outputshapefile=getfieldvalue(options,'outputshapefile','');
    175172
    176                         if (exist(options,'epsgshapefile') &  exist(options,'projshapefile')),
     173                        if (exist(options,'epsgshapefile') & exist(options,'projshapefile')),
    177174                                error('Basin shapefilecrop error message: cannot specify epsgshapefile and projshapefile at the same time');
    178175                        end
    179176                        if exist(options,'epsgshapefile'),
    … …  
    185182                        end
    186183
    187184
    188                         %create list of contours that have critical length > threshold:  (in lat,long)
     185                        %create list of contours that have critical length > threshold (in lat,long):
    189186                        contours=shpread(inshapefile);
    190187                        llist=[];
    191188                        for i=1:length(contours),
    … …  
    211208                        flags=zeros(length(contours),1);
    212209                        for i=1:length(contours),
    213210                                h=contours(i);
    214                                 in=inpolygon(h.x,h.y,ba.x,ba.y);
    215                                 if ~isempty(find(in==0)),
     211                                isin=inpolygon(h.x,h.y,ba.x,ba.y);
     212                                if ~isempty(find(isin==0)),
    216213                                        flags(i)=1;
    217214                                end
    218215                        end

  • ../trunk-jpl/src/m/qmu/helpers.py

     
    1 import numpy as np
    21from collections import OrderedDict
    32from copy import deepcopy
    43
     4import numpy as np
    55
     6
    67class struct(object):
    78    '''An empty struct that can be assigned arbitrary attributes'''
    89    pass
    … …  
    1011
    1112class Lstruct(list):
    1213    '''
    13 An empty struct that can be assigned arbitrary attributes
    14     but can also be accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
     14    An empty struct that can be assigned arbitrary attributes but can also be
     15    accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
    1516
    16 Note that 'x' returns the array and x.__dict__ will only return
    17     attributes other than the array
     17    Note that 'x' returns the array and x.__dict__ will only return attributes
     18    other than the array
    1819
    19 List-based and struct-based behaviors work normally, however they are referenced
    20     as if the other does not exist; len(x) corresponds only to
    21     the list component of x, len(x.a) corresponds to x.a, x.__dict__
    22     corresponds only to the non-x-list attributes
     20    List-based and struct-based behaviors work normally, however they are
     21    referenced as if the other does not exist; len(x) corresponds only to the
     22    list component of x, len(x.a) corresponds to x.a, x.__dict__ corresponds
     23    only to the non-x-list attributes
    2324
    24 Example uses:
     25    Examples:
    2526
    26 x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4]
    27     x.a = 'hello'
    28     len(x) -> 4
    29     x.append(5)
    30     len(x) -> 5
    31     x[2] -> 3
    32     x.a -> 'hello'
    33     print x -> [1, 2, 3, 4, 5]
    34     x.__dict__ -> {'a':'hello'}
    35     x.b = [6, 7, 8, 9]
    36     x.b[-1] -> 9
    37     len(x.b) -> 4
     27        x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4]
     28        x.a = 'hello'
     29        len(x) -> 4
     30        x.append(5)
     31        len(x) -> 5
     32        x[2] -> 3
     33        x.a -> 'hello'
     34        print x -> [1, 2, 3, 4, 5]
     35        x.__dict__ -> {'a':'hello'}
     36        x.b = [6, 7, 8, 9]
     37        x.b[-1] -> 9
     38        len(x.b) -> 4
    3839
    39 Other valid constructors:
    40           x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3]
    41           x = Lstruct(1, 2, 3)(a = 'hello')
    42           x = Lstruct([1, 2, 3], x = 'hello')
    43           x = Lstruct((1, 2, 3), a = 'hello')
     40    Other valid constructors:
     41        x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3]
     42        x = Lstruct(1, 2, 3)(a = 'hello')
     43        x = Lstruct([1, 2, 3], x = 'hello')
     44        x = Lstruct((1, 2, 3), a = 'hello')
    4445
    45 Credit: https://github.com/Vectorized/Python-Attribute-List
    46 '''
     46    Sources:
     47    -https://github.com/Vectorized/Python-Attribute-List
     48    '''
    4749
    4850    def __new__(self, *args, **kwargs):
    4951        return super(Lstruct, self).__new__(self, args, kwargs)
    … …  
    6264
    6365class OrderedStruct(object):
    6466    '''
    65 A form of dictionary-like structure that maintains the
    66     ordering in which its fields/attributes and their
    67     corresponding values were added.
     67    A form of dictionary-like structure that maintains the ordering in which
     68    its fields/attributes and their corresponding values were added.
    6869
    69 OrderedDict is a similar device, however this class
    70     can be used as an "ordered struct/class" giving
    71     it much more flexibility in practice. It is
     70    OrderedDict is a similar device, however this class can be used as an
     71    "ordered struct/class" giving it much more flexibility in practice. It is
    7272    also easier to work with fixed valued keys in-code.
    7373
    74 Eg:
    75 OrderedDict:  # a bit clumsy to use and look at
    76     x['y'] = 5
     74    Example:
     75        OrderedDict:  # a bit clumsy to use and look at
     76            x['y'] = 5
    7777
    78 OrderedStruct:  # nicer to look at, and works the same way
    79     x.y = 5
    80     OR
    81     x['y'] = 5  # supports OrderedDict-style usage
     78        OrderedStruct:  # nicer to look at, and works the same way
     79            x.y = 5
     80            OR
     81            x['y'] = 5  # supports OrderedDict-style usage
    8282
    83 Supports: len(x), str(x), for-loop iteration.
    84 Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
     83    Supports: len(x), str(x), for-loop iteration.
     84    Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
    8585
    86 Usage:
    87     x = OrderedStruct()
    88     x.y = 5
    89     x.z = 6
    90     OR
    91     x = OrderedStruct('y', 5, 'z', 6)
     86    Usage:
     87        x = OrderedStruct()
     88        x.y = 5
     89        x.z = 6
     90        OR
     91        x = OrderedStruct('y', 5, 'z', 6)
    9292
    93     # note below that the output fields as iterables are always
    94     #    in the same order as the inputs
     93    note below that the output fields as iterables are always in the same
     94    order as the inputs
    9595
    9696    x.keys() -> ['y', 'z']
    9797    x.values() -> [5, 6]
    … …  
    110110    ('x', 5)
    111111    ('y', 6)
    112112
    113     Note: to access internal fields use dir(x)
    114         (input fields will be included, but
    115         are not technically internals)
     113    Note: to access internal fields use dir(x) (input fields will be included,
     114    but are not technically internals)
    116115    '''
    117116
    118117    def __init__(self, *args):
    119         '''Provided either nothing or a series of strings, construct a structure that will,
    120 when accessed as a list, return its fields in the same order in which they were provided'''
     118        '''
     119        Provided either nothing or a series of strings, construct a structure
     120        that will, when accessed as a list, return its fields in the same order
     121        in which they were provided
     122        '''
    121123
    122     # keys and values
     124        # keys and values
    123125        self._k = []
    124126        self._v = []
    125127
    … …  
    184186            yield(a, b)
    185187
    186188    def __copy__(self):
    187         # shallow copy, hard copies of trivial attributes,
    188         # references to structures like lists/OrderedDicts
    189         # unless redefined as an entirely different structure
     189        '''
     190        shallow copy, hard copies of trivial attributes,
     191        references to structures like lists/OrderedDicts
     192        unless redefined as an entirely different structure
     193        '''
    190194        newInstance = type(self)()
    191195        for k, v in list(self.items()):
    192196            exec(('newInstance.%s = v') % (k))
    … …  
    193197        return newInstance
    194198
    195199    def __deepcopy__(self, memo=None):
    196         # hard copy of all attributes
    197         # same thing but call deepcopy recursively
    198         # technically not how it should be done,
    199         # (see https://docs.python.org/2/library/copy.html  #copy.deepcopy )
    200         # but will generally work in this case
     200        '''
     201        hard copy of all attributes
     202        same thing but call deepcopy recursively
     203        technically not how it should be done,
     204        (see https://docs.python.org/2/library/copy.html  #copy.deepcopy )
     205        but will generally work in this case
     206        '''
    201207        newInstance = type(self)()
    202208        for k, v in list(self.items()):
    203209            exec(('newInstance.%s = deepcopy(v)') % (k))
    … …  
    211217        i = self._k.index(key)
    212218        k = self._k.pop(i)
    213219        v = self._v.pop(i)
    214     #exec('del self.%s')%(key)
     220        #exec('del self.%s')%(key)
    215221        return (k, v)
    216222
    217223    def keys(self):
    … …  
    226232
    227233def isempty(x):
    228234    '''
    229     returns true if object is +  -infinity, NaN, None, '', has length 0, or is an
    230     array/matrix composed only of such components (includes mixtures of "empty" types)'''
     235    returns true if object is +/-infinity, NaN, None, '', has length 0, or is
     236    an array/matrix composed only of such components (includes mixtures of
     237    "empty" types)
     238    '''
    231239
    232240    if type(x) in [list, np.ndarray, tuple]:
    233241        if np.size(x) == 0:
    … …  
    261269
    262270
    263271def fieldnames(x, ignore_internals=True):
    264     '''returns a list of fields of x
    265     ignore_internals ignores all fieldnames starting with '_' and is True by default'''
     272    '''
     273    returns a list of fields of x
     274    ignore_internals ignores all fieldnames starting with '_' and is True by
     275    default
     276    '''
    266277    result = list(vars(x).keys())
    267278
    268279    if ignore_internals:
    … …  
    272283
    273284
    274285def isfield(x, y, ignore_internals=True):
    275     '''is y is a field of x?
    276     ignore_internals ignores all fieldnames starting with '_' and is True by default'''
     286    '''
     287    is y is a field of x?
     288    ignore_internals ignores all fieldnames starting with '_' and is True by
     289    default
     290    '''
    277291    return str(y) in fieldnames(x, ignore_internals)
    278292
    279293
    … …  
    280294def fileparts(x):
    281295    '''
    282296    given:   "path/path/.../file_name.ext"
    283     returns: [path, file_name, ext] (list of strings)'''
     297    returns: [path, file_name, ext] (list of strings)
     298    '''
    284299    try:
    285300        a = x[:x.rindex('/')]  #path
    286301        b = x[x.rindex('/') + 1:]  #full filename
    … …  
    296311
    297312def fullfile(*args):
    298313    '''
    299     use:
     314    usage:
     315        fullfile(path, path, ... , file_name + ext)
    300316
    301     fullfile(path, path, ... , file_name + ext)
    302317    returns: "path/path/.../file_name.ext"
    303318
    304319    with all arguments as strings with no "/"s
    305320
    306321    regarding extensions and the '.':
    307     as final arguments ('file.doc') or ('file' + '.doc') will work
    308     ('final', '.doc'), and the like, will not (you'd get 'final/.doc')
    309 '''
     322        as final arguments ('file.doc') or ('file' + '.doc') will work
     323        ('final', '.doc'), and the like, will not (you'd get 'final/.doc')
     324    '''
    310325    result = str(args[0])
    311326    for i in range(len(args[1:])):
    312327        # if last argument wasn't empty, add a '/' between it and the next argument
    … …  
    320335def findline(fidi, s):
    321336    '''
    322337    returns full first line containing s (as a string), or None
     338
    323339    Note: will include any newlines or tabs that occur in that line,
    324     use str(findline(f, s)).strip() to remove these, str() in case result is None'''
     340    use str(findline(f, s)).strip() to remove these, str() in case result is
     341    None
     342    '''
    325343    for line in fidi:
    326344        if s in line:
    327345            return line
    … …  
    330348
    331349def empty_nd_list(shape, filler=0., as_numpy_ndarray=False):
    332350    '''
    333 returns a python list of the size/shape given (shape must be int or tuple)
     351    returns a python list of the size/shape given (shape must be int or tuple)
    334352    the list will be filled with the optional second argument
    335353
    336354    filler is 0.0 by default
    337355
    338     as_numpy_ndarray will return the result as a numpy.ndarray and is False by default
     356    as_numpy_ndarray will return the result as a numpy.ndarray and is False by
     357    default
    339358
    340     Note: the filler must be either None/np.nan/float('NaN'), float/double, or int
    341         other numpy and float values such as +/- np.inf will also work
     359    Note: the filler must be either None/np.nan/float('NaN'), float/double, or
     360    int. other numpy and float values such as +/- np.inf will also work
    342361
    343 use:
    344     empty_nd_list((5, 5), 0.0)  # returns a 5x5 matrix of 0.0's
    345     empty_nd_list(5, None)  # returns a 5 long array of NaN
    346 '''
     362    Usage:
     363        empty_nd_list((5, 5), 0.0)  # returns a 5x5 matrix of 0.0's
     364        empty_nd_list(5, None)  # returns a 5 long array of NaN
     365    '''
    347366    result = np.empty(shape)
    348367    result.fill(filler)
    349368    if not as_numpy_ndarray:

  • ../trunk-jpl/src/m/mech/analyticaldamage.py

     
    11import numpy as np
     2
    23from averaging import averaging
    3 #from plotmodel import plotmodel
    44from thomasparams import thomasparams
    55
    66

  • ../trunk-jpl/src/m/plot/applyoptions.m

     
    307307        end
    308308end
    309309
    310 
    311310%text (default value is empty, not NaN...)
    312311if exist(options,'text');
    313312        textstring=getfieldvalue(options,'text');
    … …  
    351350        scaleruler(options);
    352351end
    353352
    354 %streamliness
     353%streamlines
    355354if exist(options,'streamlines'),
    356355        plot_streamlines(md,options);
    357356end

  • ../trunk-jpl/src/m/plot/plot_manager.py

     
    1 
    2 from checkplotoptions import checkplotoptions
    3 from plot_mesh import plot_mesh
    4 from plot_BC import plot_BC
    5 from plot_elementnumbering import plot_elementnumbering
    6 from plot_vertexnumbering import plot_vertexnumbering
    7 from processmesh import processmesh
    8 from processdata import processdata
    9 from plot_unit import plot_unit
    10 from applyoptions import applyoptions
    11 
    121try:
    132    from osgeo import gdal
    143    overlaysupport = True
    … …  
    165    print('osgeo/gdal for python not installed, overlay plots are not enabled')
    176    overlaysupport = False
    187
     8from applyoptions import applyoptions
     9from checkplotoptions import checkplotoptions
     10from plot_BC import plot_BC
     11from plot_mesh import plot_mesh
     12from plot_elementnumbering import plot_elementnumbering
    1913if overlaysupport:
    2014    from plot_overlay import plot_overlay
     15from plot_unit import plot_unit
     16from plot_vertexnumbering import plot_vertexnumbering
     17from processdata import processdata
     18from processmesh import processmesh
    2119
    2220
    2321def plot_manager(md, options, fig, axgrid, gridindex):
    … …  
    2624
    2725    'fig' is a handle to the figure instance created by plotmodel.
    2826
    29     'ax' is a handle to the axes instance created by plotmodel. This is
     27    'axgrid' is a handle to the axes instance created by plotmodel. This is
    3028    currently generated using matplotlib's AxesGrid toolkit.
    3129
     30    'gridindex' is passed through to specialized plot* functions and
     31    applyoptions.
     32
    3233    Usage:
    3334        plot_manager(md, options, fig, ax)
    3435
    … …  
    3637    '''
    3738    #parse options and get a structure of options
    3839    options = checkplotoptions(md, options)
     40
    3941    #get data to be plotted
    4042    data = options.getfieldvalue('data')
     43
    4144    #add ticklabel has a default option
     45    #
     46    # TODO: Check why we are doing this and if it is absolutely necessary (see
     47    #       also src/m/plot/plot_mesh.py, src/m/plot/applyoptions.py)
    4248    options.addfielddefault('ticklabels', 'on')
    4349
    4450    ax = axgrid[gridindex]
    … …  
    5763    if isinstance(data, str):
    5864        if data == 'mesh':
    5965            plot_mesh(md, options, fig, axgrid, gridindex)
    60 
    61     #fig.delaxes(fig.axes[1])  # hack to remove colorbar after the fact
     66            #fig.delaxes(fig.axes[1])  # hack to remove colorbar after the fact
    6267            return
    6368        elif data == 'BC':
    6469            plot_BC(md, options, fig, axgrid, gridindex)

  • ../trunk-jpl/src/m/plot/plot_overlay.py

     
    1 import numpy as np
    2 from processmesh import processmesh
    3 from processdata import processdata
    4 from xy2ll import xy2ll
     1import os
     2
    53import matplotlib.pyplot as plt
    64import matplotlib as mpl
    7 import os
    85try:
    96    from mpl_toolkits.basemap import Basemap
    107except ImportError:
    118    print('Basemap toolkit not installed')
     9import numpy as np
    1210try:
    1311    from osgeo import gdal
    1412except ImportError:
    1513    print('osgeo/gdal for python not installed, plot_overlay is disabled')
    1614
     15from processdata import processdata
     16from processmesh import processmesh
     17from xy2ll import xy2ll
    1718
     19
    1820def plot_overlay(md, data, options, ax):
    1921    '''
    20     Function for plotting a georeferenced image.  This function is called
    21     from within the plotmodel code.
     22    Function for plotting a georeferenced image. Called from plot_manager by
     23    call to plotmodel.
     24
     25    Usage:
     26        plot_overlay(md, data, options, ax)
     27
     28    See also: PLOTMODEL
    2229    '''
    2330
    2431    x, y, z, elements, is2d, isplanet = processmesh(md, [], options)
    … …  
    8087        plt.hist(arr.flatten(), bins=256, range=(0., 1.))
    8188        plt.title('histogram of overlay image, use for setting overlaylims')
    8289        plt.show()
    83         plt.sca(ax)  # return to original axes / figure
     90        plt.sca(ax) # return to original axes/figure
    8491
    8592    # get parameters from cropped geotiff
    8693    trans = gtif.GetGeoTransform()

  • ../trunk-jpl/src/m/plot/plotmodel.m

     
    99numberofplots=options.numberofplots;
    1010
    1111%get number of subplots
    12 subplotwidth=ceil(sqrt(options.numberofplots));
     12subplotwidth=ceil(sqrt(numberofplots));
    1313
    1414%if nlines and ncols specified, then bypass.
    1515if ~exist(options.list{1},'nlines') & ~exist(options.list{1},'ncols')
    … …  
    103103                end
    104104        end % }}}
    105105
    106         %Use zbuffer renderer (snoother colors) and white background
     106        %Use zbuffer renderer (smoother colors) and white background
    107107        set(f,'Renderer','zbuffer','color',getfieldvalue(options.list{1},'figurebackgroundcolor','w'));
    108108
    109109        %Go through all data plottable and close window if an error occurs
    … …  
    124124                rethrow(me);
    125125        end
    126126else
    127         error('plotmodel error message: no output data found. ');
     127        error('plotmodel error message: no output data found.');
    128128end

  • ../trunk-jpl/src/m/plot/plot_manager.m

     
    1414
    1515%figure out if this is a special plot
    1616if ischar(data),
    17 
    1817        switch data,
    19 
    2018                case 'boundaries',
    2119                        plot_boundaries(md,options,subplotwidth,i);
    2220                        return;
    … …  
    3230                case 'highlightelements',
    3331                        plot_highlightelements(md,options,subplotwidth,i);
    3432                        return;
    35 
    3633                case 'qmumean',
    3734                        plot_qmumean(md,options,nlines,ncols,i);
    3835                        return;
    39 
    4036                case 'qmustddev',
    4137                        plot_qmustddev(md,options,nlines,ncols,i);
    4238                        return;
    43 
    4439                case 'qmuhistnorm',
    4540                        plot_qmuhistnorm(md,options,nlines,ncols,i);
    4641                        return;
    47 
    4842                case 'qmu_mass_flux_segments',
    4943                        plot_qmu_mass_flux_segments(md,options,nlines,ncols,i);
    5044                        return;
    51 
    5245                case 'part_hist',
    5346                        plot_parthist(md,options,nlines,ncols,i);
    5447                        return;
    … …  
    120113                case 'segments'
    121114                        plot_segments(md,options,subplotwidth,i,data)
    122115                        return
    123 
    124116                case 'quiver'
    125117                        data=[md.initialization.vx md.initialization.vy]; %Go ahead and try plot_unit
    126 
    127118                case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',...
    128119                                'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',...
    129120                                'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'},
    … …  
    137128                        return;
    138129                case 'transient_results',
    139130                        plot_transient_results(md,options,subplotwidth,i);
    140 
    141131                case 'transient_field',
    142132                        plot_transient_field(md,options,subplotwidth,i);
    143133                        return;
    144 
    145134        otherwise,
    146 
    147135                if ismember(data,properties('model')),
    148136                        data=eval(['md.' data ';']);
    149137                else
    … …  
    158146        return;
    159147end
    160148
    161 %Figure out if this is a semi-transparent plot.
     149%Figure out if this is a Google Maps plot.
    162150if exist(options,'googlemaps'),
    163151        plot_googlemaps(md,data,options,nlines,ncols,i);
    164152        return;
    165153end
    166154
    167 %Figure out if this is a semi-transparent plot.
     155%Figure out if this is a landsat plot.
    168156if getfieldvalue(options,'landsat',0),
    169157        plot_landsat(md,data,options,nlines,ncols,i);
    170158        return;
    171159end
    172160
    173 %Figure out if this is a semi-transparent plot.
     161%Figure out if this is a gridded plot.
    174162if exist(options,'gridded'),
    175163        plot_gridded(md,data,options,nlines,ncols,i);
    176164        return;

  • ../trunk-jpl/src/m/shp/shpread.py

     
     1from os import path
     2try:
     3    import shapefile
     4except ImportError:
     5    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
     6
     7from pairoptions import pairoptions
     8
     9
     10def shpread(filename, *args): #{{{
     11    '''
     12    SHPREAD - read a shape file and build a dict
     13
     14    This routine reads a shape file .shp and builds a dict array containing the
     15    fields x and y corresponding to the coordinates, one for the filename
     16    of the shp file, for the density, for the nodes, and a field closed to
     17    indicate if the domain is closed. If this initial shapefile is point only,
     18    the fields closed and points are ommited.
     19    The first argument is the .shp file to be read and the second one
     20    (optional) indicates if the last point shall be read (1 to read it, 0 not
     21    to).
     22
     23    The current implementation of shpread depends on PyShp.
     24
     25    Usage:
     26        dict = shpread(filename)
     27
     28    Example:
     29        From underling PyShp implementation, "The shapefile format is actually
     30        a collection of three files. You specify the base filename of the
     31        shapefile or the complete filename of any of the shapefile component
     32        files.""
     33
     34        dict = shpread('domainoutline.shp')
     35        OR
     36        dict = shpread('domainoutline.dbf')
     37        OR
     38        dict = shpread('domainoutline')
     39
     40        "OR any of the other 5+ formats which are potentially part of a
     41        shapefile. The library does not care about file extensions". We do,
     42        however, check that a file with the base filename or base filename with
     43        .shp extension exists.
     44
     45    See also EXPDOC, EXPWRITEASVERTICES
     46
     47    Sources:
     48    - https://github.com/GeospatialPython/pyshp
     49    '''
     50
     51    #recover options
     52    options = pairoptions(*args)
     53
     54    #some checks
     55    if not (path.isfile(filename) or path.isfile(filename + '.shp')):
     56        raise RuntimeError('shpread error message: file {} or {}.shp not found!'.format(filename, filename))
     57
     58    #read shapefile
     59    sf = shapefile.Reader(filename)
     60
     61    Dicts = []
     62
     63    #retrieve ID (if it exists)
     64    fields = sf.fields
     65    for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
     66        if fields[i][0] == 'id':
     67            name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
     68            break
     69
     70    shapes = sf.shapes()
     71    for i in range(len(shapes)):
     72        Dict = {}
     73        shape = shapes[i]
     74        if shape.shapeType == shapefile.POINT:
     75            Dict['x'] = shape.points[0][0]
     76            Dict['y'] = shape.points[0][1]
     77            Dict['density'] = 1
     78            Dict['Geometry'] = 'Point'
     79        elif shape.shapeType == shapefile.POLYLINE:
     80            num_points = len(shape.points)
     81            x = []
     82            y = []
     83            for j in range(num_points):
     84                point = shape.points[j]
     85                x.append(point[0])
     86                y.append(point[1])
     87            Dict['x'] = x
     88            Dict['y'] = y
     89            Dict['nods'] = num_points
     90            Dict['density'] = 1
     91            Dict['closed'] = 1
     92            Dict['BoundingBox'] = shape.bbox
     93            Dict['Geometry'] = 'Line'
     94        elif shape.shapeType == shapefile.POLYGON:
     95            num_points = len(shape.points)
     96            x = []
     97            y = []
     98            for j in range(num_points):
     99                point = shape.points[j]
     100                x.append(point[0])
     101                y.append(point[1])
     102            Dict['x'] = x
     103            Dict['y'] = y
     104            Dict['nods'] = num_points
     105            Dict['density'] = 1
     106            Dict['closed'] = 1
     107            Dict['BoundingBox'] = shape.bbox
     108            Dict['Geometry'] =  'Polygon'
     109        else:
     110            # NOTE: We could do this once before looping over shapes as all
     111            #       shapes in the file must be of the same type, but we would
     112            #       need to have a second check anyway in order to know how to
     113            #       parse the points. So, let's just assume the file is not
     114            #       malformed.
     115            #
     116            raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName))
     117        name = ''
     118        fields = sf.fields
     119        for j in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
     120            fieldname = fields[j][0]
     121            # 'id' field gets special treatment
     122            if fieldname == 'id':
     123                name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
     124            else:
     125                Dict[fieldname] = sf.record(j - 1)[0]
     126        Dict['name'] = name
     127        Dicts.append(Dict)
     128
     129    invert = options.getfieldvalue('invert', 0)
     130    if invert:
     131        for i in range(len(Dicts)):
     132            Dicts[i].x = np.flipud(Dicts[i].x)
     133            Dicts[i].y = np.flipud(Dicts[i].y)
     134
     135    return Dicts
     136#}}}

  • ../trunk-jpl/src/m/shp/shpread.m

     
    11function Struct=shpread(filename,varargin)
    2 %SHPREAD - read a shape file and build a Structure
     2%SHPREAD - read a shape file and build a struct
    33%
    4 %   This routine reads a shape file .shp and builds a Structure containing the
    5 %   fields x and y corresponding to the coordinates, one for the filename of
    6 %   the shp file, for the density, for the nodes, and a field closed to
    7 %   indicate if the domain is closed.
    8 %   If this initial shapefile is point only, the fields closed and
    9 %   points are ommited
    10 %   The first argument is the .shp file to be read and the second one (optional)
    11 %   indicates if the last point shall be read (1 to read it, 0 not to).
     4%       This routine reads a shape file .shp and builds a structure array
     5%       containing the fields x and y corresponding to the coordinates, one for the
     6%       filename of the shp file, for the density, for the nodes, and a field
     7%       closed to indicate if the domain is closed. If this initial shapefile is
     8%       point only, the fields closed and points are omitted.
     9%       The first argument is the .shp file to be read and the second one
     10%       (optional) indicates if the last point shall be read (1 to read it, 0 not
     11%       to).
    1212%
    13 %   Usage:
    14 %      Struct=shpread(filename)
     13%       Usage:
     14%               Struct=shpread(filename)
    1515%
    16 %   Example:
    17 %      Struct=shpread('domainoutline.shp')
     16%       Example:
     17%               Struct=shpread('domainoutline.shp')
    1818%
    19 %   See also EXPDOC, EXPWRITEASVERTICES
     19%       See also EXPDOC, EXPWRITEASVERTICES
    2020
    2121%recover options
    2222options=pairoptions(varargin{:});
    … …  
    2626        error(['shpread error message: file ' filename ' not found!']);
    2727end
    2828
    29 %initialize number of profile
    30 count=0;
    31 
    3229%read shapefile
    3330shp=shaperead(filename);
    3431
    … …  
    3532Struct=struct([]);
    3633fields=fieldnames(shp);
    3734for i=1:length(shp),
    38         if strcmpi(shp(i).Geometry,'Polygon'),
     35        if strcmpi(shp(i).Geometry,'Point'),
    3936                x=shp(i).X'; y=shp(i).Y';
    4037                ids=find(isnan(x));
    4138                x(ids)=[]; y(ids)=[];
    … …  
    4239
    4340                Struct(end+1).x=x;
    4441                Struct(end).y=y;
    45                 Struct(end).nods=length(x);
    4642                Struct(end).density=1;
    47                 Struct(end).closed=1;
    4843                if isfield(shp,'id'),
    4944                        Struct(end).name=num2str(shp(i).id);
    5045                else
    … …  
    5651                                Struct(end).(field)=shp(i).(field);
    5752                        end
    5853                end
    59         end
    60        
    61         if strcmpi(shp(i).Geometry,'Line'),
     54        elseif strcmpi(shp(i).Geometry,'Line'),
    6255                x=shp(i).X'; y=shp(i).Y';
    6356                ids=find(isnan(x));
    6457                x(ids)=[]; y(ids)=[];
    … …  
    7972                                Struct(end).(field)=shp(i).(field);
    8073                        end
    8174                end
    82         end
    83 
    84 
    85         if strcmpi(shp(i).Geometry,'Point'),
     75        elseif strcmpi(shp(i).Geometry,'Polygon'),
    8676                x=shp(i).X'; y=shp(i).Y';
    8777                ids=find(isnan(x));
    8878                x(ids)=[]; y(ids)=[];
    … …  
    8979
    9080                Struct(end+1).x=x;
    9181                Struct(end).y=y;
     82                Struct(end).nods=length(x);
    9283                Struct(end).density=1;
     84                Struct(end).closed=1;
    9385                if isfield(shp,'id'),
    9486                        Struct(end).name=num2str(shp(i).id);
    9587                else

  • ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py

     
    1 import os
    21import subprocess
    32
    43import numpy as np
    … …  
    98from pairoptions import *
    109
    1110
    12 def gmshplanet(*varargin):
     11def gmshplanet(*args):
    1312    '''
    1413    GMSHPLANET - mesh generation for a sphere. Very specific code for gmsh from $ISSM_DIR/src/demos/simple_geo/sphere.geo
    1514
    … …  
    3332        raise RuntimeError('gmshplanet: Gmsh major version %s not supported!' % gmshmajorversion)
    3433
    3534    #process options
    36     options = pairoptions(*varargin)
     35    options = pairoptions(*args)
    3736    #options = deleteduplicates(options, 1)
    3837
    3938    #recover parameters:

  • ../trunk-jpl/src/m/coordsystems/gdaltransform.py

     
    55
    66from loadvars import *
    77
     8
    89def gdaltransform(x, y, proj_in, proj_out): #{{{
    910    '''
    1011    GDALTRANSFORM - switch from one projection system to another
    1112
    12         Usage:
    13             [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out)
     13    Usage:
     14        [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out)
    1415
    15         Example:
    16             [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031')
     16    Example:
     17        [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031')
    1718
    18         For reference:
    19             EPSG: 4326 (lat, long)
    20             EPSG: 3341 (Greenland,  UPS 45W, 70N)
    21             EPSG: 3031 (Antarctica, UPS 0E,  71S)
     19    For reference:
     20        EPSG: 4326 (lat, long)
     21        EPSG: 3341 (Greenland,  UPS 45W, 70N)
     22        EPSG: 3031 (Antarctica, UPS 0E,  71S)
    2223
    23         ll2xy default projection Antarctica:
    24             +proj = stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
    25         ll2xy default projection Greenland:
    26             +proj = stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
    27         Bamber's Greenland projection
    28             +proj = stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     24    ll2xy default projection Antarctica:
     25        +proj = stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     26    ll2xy default projection Greenland:
     27        +proj = stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
     28    Bamber's Greenland projection
     29        +proj = stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
    2930
    30         To get proj.4 string form EPSG, use gdalsrsinfo. Example:
     31    To get proj.4 string form EPSG, use gdalsrsinfo. Example:
    3132
    32             gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
     33        gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
     34
     35    TODO:
     36    - Remove shlex and follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py
    3337    '''
    3438
    3539    # Give ourselves unique file names
    … …  
    4145    file_in.write(b'%8g %8g\n' % (x.flatten(1) y.flatten(1)))
    4246    file_in.close()
    4347
    44     args = shlex.split('gdaltransform -s_srs %s -t_srs %s < %s > %s' % (proj_in, proj_out, filename_in, filename_out))
     48    args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out))
    4549    subprocess.check_call(args) # Will raise CalledProcessError if return code is not 0
    4650
    4751    A = loadvars(filename_out)
    … …  
    5256
    5357    os.remove(filename_in)
    5458    os.remove(filename_out)
     59
     60    return [xout, yout]
    5561#}}}

  • ../trunk-jpl/src/m/geometry/polyarea.py

     
     1import math
     2
     3import numpy as np
     4
     5
     6def polyarea: #{{{
     7    '''
     8    POLYAREA - returns the area of the 2-D polygon defined by the vertices in
     9    lists x and y
     10   
     11    Partial implementation of MATLAB's polyarea function. If x and y are
     12    lists of the same length, then polyarea returns the scalar area of the
     13    polygon defined by x and y.
     14
     15    Usage:
     16        a = polyarea(x, y)
     17
     18    Sources:
     19    - https://www.mathworks.com/help/matlab/ref/polyarea.html
     20    - https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates
     21    - https://en.wikipedia.org/wiki/Shoelace_formula
     22
     23    TODO:
     24    - Test that output falls within some tolerance of MATLAB's polyarea
     25    function.
     26    '''
     27    return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
     28#}}}

  • ../trunk-jpl/src/m/classes/boundary.py

     
     1import math
     2
     3import numpy as np
     4
     5from epsg2proj import epsg2proj
     6from pairoptions import pairoptions
     7from shpread import shpread
     8
     9
     10class boundary(object): #{{{
     11    '''
     12    BOUNDARY class definition
     13
     14    Usage:
     15        boundary = boundary()
     16    '''
     17
     18    def __init__(self, *args): #{{{
     19        self.boundaries = []
     20        self.name       = ''
     21        self.continent  = ''
     22        self.proj       = epsg2proj(4326)
     23
     24        self.setdefaultparameters()
     25
     26        if len(args):
     27            options = pairoptions(*args)
     28
     29            #recover field values
     30            self.shppath = options.getfieldvalue('shppath', '')
     31            self.shpfilename = options.getfieldvalue('shpfilename', '')
     32            self.orientation = options.getfieldvalue('orientation', 'normal')
     33
     34            #figure out projection string:
     35            if options.exist('epsg'):
     36                if options.exist('proj'):
     37                    raise RuntimeError('Error in constructor for boundary {}. Cannot supply epsg and proj at the same time!'.format(self.shppath))
     38                epsg = options.getfieldvalue('epsg')
     39                proj = epsg2proj(epsg) # convert to PROJ.4 format
     40            elif options.exist('proj'):
     41                if options.exist('epsg'):
     42                    raise RuntimeError('Error in constructor for boundary {}. Cannot supply proj and epsg at the same time!'.format(self.shppath))
     43                proj = options.getfieldvalue('proj')
     44            else:
     45                proj = '+proj=longlat +datum=WGS84 +no_defs'
     46
     47            self.proj = proj
     48    #}}}
     49
     50    def __repr__(self): #{{{
     51        s = '   boundary parameters:\n'
     52        s += '{}\n'.format(fielddisplay(self, 'shppath', 'path to filename for this boundary'))
     53        s += '{}\n'.format(fielddisplay(self, 'shpfilename', 'shape file name'))
     54        s += '{}\n'.format(fielddisplay(self, 'orientation', 'orientation (default is \'normal\', can be \'reverse\')'))
     55        s += '{}\n'.format(fielddisplay(self, 'proj', 'shape file projection string (follows PROJ.4 standard)'))
     56
     57        return s
     58    #}}}
     59
     60    def setdefaultparameters(self): #{{{
     61        self.shppath = ''
     62        self.shpfilename = ''
     63        self.orientation = 'normal'
     64        self.proj = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326
     65
     66        return self
     67    #}}}
     68
     69    def name(self): #{{{
     70        output = self.shpfilename
     71
     72        return output
     73    #}}}
     74
     75    def edges(self): #{{{
     76        #read domain
     77        path, name, ext = fileparts(shp.shpfilename)
     78        if ext != 'shp':
     79            ext = 'shp'
     80        output = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     81
     82        #do we reverse?
     83        if self.orientation == 'reverse':
     84            output.x = np.flipud(output.x)
     85            output.y = np.flipud(output.y)
     86    #}}}
     87
     88    def plot(self, *args): #{{{
     89        #recover options
     90        options = pairoptions(*args)
     91
     92        #parse input
     93        # TODO: Sort out what options we should set (see
     94        # src/m/classes/boundary.m)
     95       
     96        if options.exist('epsg'):
     97            if options.exist('proj'):
     98                raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection')
     99            proj = epsg2proj(options.getfieldvalue('epsg'))
     100        elif options.exist('proj'):
     101            proj = options.getfieldvalue('proj')
     102        else:
     103            proj = epsg2proj(4326)
     104
     105        #read domain
     106        path, name, ext = fileparts(shp.shpfilename)
     107        if ext != 'shp':
     108            ext = 'shp'
     109        domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     110
     111        #convert boundary to another reference frame #{{{
     112        for i in range(len(domain)):
     113            try:
     114                x, y = gdaltransform(domain[i].x, domain[i].y, self.proj, proj)
     115            except error as e:
     116                print(e)
     117                print(self)
     118
     119        # TODO: Figure out how to recover figure here: do we pass 'fig' and
     120        # 'ax' in args?
     121        #for i in range(len(domain)):
     122            # x = domain[i].x * unitmultiplier
     123            # y = domain[i].y * unitmultiplier
     124            # if len(x) == 1:
     125        #}}}
     126
     127        #TODO: Finish translating from MATLAB after test2004.py runs without plot
     128    #}}}
     129
     130    def checkconsistency(self, *args): #{{{
     131        #recover options
     132        options = pairoptions(*args)
     133        tolerance = getfieldvalue(options, 'tolerance', 1e-5)
     134
     135        #recover edges
     136        edges = self.edges()
     137
     138        if edges.Geometry == 'Point':
     139            return
     140        else:
     141            #check that we don't have two identical vertices
     142            x = edges.x
     143            y = edges.y
     144            distance = math.sqrt((x[:-2] - x[1:-1]) ** 2 + (y[:-2] - y[1:-1]) ** 2)
     145            dmax = distance.max()
     146            isidentical = np.where(np.asarray(distance) < dmax * tolerance)
     147            for elem in isidentical: # distance should be a 1D array, but return from np.where is a tuple of arrays
     148                if len(elem) != 0:
     149                    raise Exception('boundary {} has two vertices extermely close to one another'.format(shpfilename))
     150
     151    def plot3d(self, *args): #{{{
     152        #recover options
     153        options = pairoptions(*args)
     154
     155        #parse input
     156        # TODO: Sort out what options we should set (see
     157        # src/m/classes/boundary.m)
     158       
     159        if options.exist('epsg'):
     160            if options.exist('proj'):
     161                raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection')
     162            proj = epsg2proj(options.getfieldvalue('epsg'))
     163        elif options.exist('proj'):
     164            proj = options.getfieldvalue('proj')
     165        else:
     166            proj = epsg2proj(4326)
     167
     168        #read domain
     169        path, name, ext = fileparts(shp.shpfilename)
     170        if ext != 'shp':
     171            ext = 'shp'
     172        domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
     173
     174        #TODO: Finish translating from MATLAB after test2004.py runs without plot
     175    #}}}
     176#}}}

  • ../trunk-jpl/src/m/classes/sealevelmodel.m

     
    162162                        list=unique(list);
    163163                end % }}}
    164164                function addbasin(self,bas) % {{{
    165                 if ~strcmpi(class(bas),'basin')
    166                         error('addbasin method only takes a ''basin'' class object as input');
    167                 end;
    168                 self.basins{end+1}=bas;
     165                        if ~strcmpi(class(bas),'basin')
     166                                error('addbasin method only takes a ''basin'' class object as input');
     167                        end;
     168                        self.basins{end+1}=bas;
    169169                end % }}}
    170170                function intersections(self,varargin) % {{{
    171171

  • ../trunk-jpl/src/m/classes/sealevelmodel.py

     
    1010from private import private
    1111from TwoDToThreeD import TwoDToThreeD
    1212
     13
    1314class sealevelmodel(object):
    1415    '''
    1516    SEALEVELMODEL class definition
    1617
    17         Usage:
    18             slm = sealevelmodel(*args)
     18    Usage:
     19        slm = sealevelmodel(*args)
    1920
    20             where args is a variable list of options
     21        where args is a variable list of options
    2122
    22         Example:
    23             slm = sealevelmodel(
    24                 'icecap', md_greenland,
    25                 'icecap', md_antarctica,
    26                 'earth', md_earth
    27             )
     23    Example:
     24        slm = sealevelmodel(
     25            'icecap', md_greenland,
     26            'icecap', md_antarctica,
     27            'earth', md_earth
     28        )
    2829    '''
    2930
    3031    def __init__(self, *args): #{{{
    … …  
    5859    #}}}
    5960
    6061    def __repr__(self): # {{{
    61         s = '%s\n' % fielddisplay(self, 'icecaps', 'ice caps')
    62         s += '%s\n' % fielddisplay(self, 'earth', 'earth')
    63         s += '%s\n' % fielddisplay(self, 'settings', 'settings properties')
    64         s += '%s\n' % fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...')
    65         s += '%s\n' % fielddisplay(self, 'miscellaneous', 'miscellaneous fields')
     62        s = '{}\n'.format(fielddisplay(self, 'icecaps', 'ice caps'))
     63        s += '{}\n'.format(fielddisplay(self, 'earth', 'earth'))
     64        s += '{}\n'.format(fielddisplay(self, 'settings', 'settings properties'))
     65        s += '{}\n'.format(fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...'))
     66        s += '{}\n'.format(fielddisplay(self, 'miscellaneous', 'miscellaneous fields'))
    6667    #}}}
    6768
    6869    def setdefaultparameters(self): # {{{
    … …  
    8485        # Is the coupler turned on?
    8586        for i in range(len(slm.icecaps)):
    8687            if slm.icecaps[i].transient.iscoupler == 0:
    87                 warnings.warn('sealevelmodel checkconsistency error: icecap model %s should have the transient coupler option turned on!' % slm.icecaps[i].miscellaneous.name)
     88                warnings.warn('sealevelmodel checkconsistency error: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name))
    8889
    8990        if slm.earth.transient.iscoupler == 0:
    9091            warnings.warn('sealevelmodel checkconsistency error: earth model should have the transient coupler option turned on!')
    … …  
    9293        # Check that the transition vectors have the right size
    9394        for i in range(len(slm.icecaps)):
    9495            if slm.icecaps[i].mesh.numberofvertices != len(slm.earth.slr.transitions[i]):
    95                 raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: %i name: %s' % (i, slm.icecaps[i].miscellaneous.name))
     96                raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name))
    9697
    9798        # Check that run frequency is the same everywhere
    9899        for i in range(len(slm.icecaps)):
    99100            if slm.icecaps[i].slr.geodetic_run_frequency != slm.earth.geodetic_run_frequency:
    100                 raise RuntimeError('sealevelmodel checkconsistency error: icecap model %s should have the same run frequency as earth!' % slm.icecaps[i].miscellaneous.name)
     101                raise RuntimeError('sealevelmodel checkconsistency error: icecap model {} should have the same run frequency as earth!'.format(slm.icecaps[i].miscellaneous.name))
    101102
    102103        # Make sure steric_rate is the same everywhere
    103104        for i in range(len(slm.icecaps)):
    104105            md = slm.icecaps[i]
    105106            if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
    106                 raise RuntimeError('steric rate on ice cap %s is not the same as for the earth' % md.miscellaneous.name)
     107                raise RuntimeError('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
    107108    #}}}
    108109
    109110    def mergeresults(self): # {{{
    … …  
    137138   
    138139    def listcaps(self): # {{{
    139140        for i in range(len(self.icecaps)):
    140             print('%i: %s' % (i, self.icecaps[i].miscellaneous.name))
     141            print('{}: {}'.format(i, self.icecaps[i].miscellaneous.name))
    141142    #}}}
    142143
    143144    def continents(self): # {{{
    … …  
    184185            yei = mdi.mesh.y[mdi.mesh.elements] * onesmatrix / 3
    185186            zei = mdi.mesh.z[mdi.mesh.elements] * onesmatrix / 3
    186187
    187             print('Computing vertex intersections for basin %s' % self.basins[i].name)
     188            print('Computing vertex intersections for basin {}'.format(self.basins[i].name))
    188189
    189190            self.transitions.append(meshintersect3d(self.earth.mesh.x, self.earth.mesh.y, self.earth.mesh.z, mdi.mesh.x, mdi.mesh.y, mdi.mesh.z, 'force', force))
    190191            self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force))

  • ../trunk-jpl/src/m/classes/basin.py

     
     1import math
     2
     3from boundary import boundary
     4from epsg2proj import epsg2proj
     5from helpers import fileparts
     6from laea import laea
     7from pairoptions import pairoptions
     8from polyarea import polyarea
     9from shpread import shpread
     10
     11
     12class basin(object): #{{{
     13    '''
     14    BASIN class definition
     15
     16    Usage:
     17        basin = basin()
     18    '''
     19
     20    def __init__(self, *args): #{{{
     21        self.boundaries = []
     22        self.name       = ''
     23        self.continent  = ''
     24        self.proj       = epsg2proj(4326)
     25
     26        self.setdefaultparameters()
     27
     28        if len(args):
     29            options = pairoptions(*args)
     30
     31            #recover field values
     32            self.boundaries = options.getfieldvalue('boundaries', [])
     33            self.name = options.getfieldvalue('name', '')
     34            self.continent = options.getfieldvalue('continent', '')
     35
     36            # figure out projection string
     37            if options.exist('epsg'):
     38                if options.exist('proj'):
     39                    raise RuntimeError('Error in constructor for basin %s. Cannot supply epsg and proj at the same time!'.format(self.name))
     40                epsg = options.getfieldvalue('epsg', '')
     41                proj = epsg2proj(epsg)#convert to PROJ.4 format
     42            elif options.exist('proj'):
     43                if options.exist('epsg'):
     44                    raise RuntimeError('Error in constructor for basin {}. Cannot supply proj and epsg at the same time!'.format(self.name))
     45                proj = options.getfieldvalue('proj', '')
     46            else:
     47                proj = '+proj=longlat +datum=WGS84 +no_defs'
     48
     49            self.proj = proj
     50    #}}}
     51
     52    def __repr__(self): # {{{
     53        s = '   basin parameters:\n'
     54        s += '{}\n'.format(fielddisplay(self, 'continent', 'continent name'))
     55        s += '{}\n'.format(fielddisplay(self, 'name', 'basin name'))
     56        s += '{}\n'.format(fielddisplay(self, 'proj', 'basin projection string (follows PROJ.4 standard'))
     57        s += '{}\n'.format(fielddisplay(self, 'boundaries','list of boundary objects'))
     58        for i in range(len(self.boundaries)):
     59            s += '             boundary #{}: {}\n'.format((i + 1), self.boundaries[i])
     60
     61        return s
     62    #}}}
     63
     64    def setdefaultparameters(self): # {{{
     65        self.name       = ''
     66        self.continent  = ''
     67        self.proj       = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326
     68        self.boundaries = []
     69
     70        return self
     71    #}}}
     72
     73    def isnameany(self, *args): #{{{
     74        for arg in args:
     75            if arg == self.name:
     76                return 1
     77        return 0
     78    #}}}
     79
     80    def iscontinentany(self, *args): #{{{
     81        for arg in args:
     82            if arg == self.continent:
     83                return 1
     84        return 0
     85    #}}}
     86
     87    def outputname(self, *args): #{{{
     88        #recover options
     89        options = pairoptions(*args)
     90        extension = options.getfieldvalue('extension', 1)
     91
     92        path, name, ext = fileparts(*args)
     93        if extension:
     94            output = '{}{}'.format(name, ext)
     95        else:
     96            output = name
     97
     98        return output
     99    #}}}
     100
     101    def plot(self, *args): #{{{
     102        #add option
     103        for i in range(len(self.boundaries)):
     104            self.boundaries[i].plot('proj', self.proj, *args)
     105    #}}}
     106
     107    def plot3d(self, *args): #{{{
     108        #add option
     109        for i in range(len(self.boundaries)):
     110            self.boundaries[i].plot3d(*args)
     111    #}}}
     112
     113    def contour(self, *args): #{{{
     114        #recover options
     115        options = pairoptions(*args)
     116        x = []
     117        y = []
     118
     119        #go through boundaries, recover edges and project them in the basin epsg reference frame
     120        for i in range(len(self.boundaries)):
     121            boundary = self.boundaries[i]
     122            contour = boundary.edges()
     123            contour.x, contour.y = gdaltransform(contour.x, contour.y, boundary.proj, self.proj)
     124            x = [[x], [contour.x]]
     125            y = [[y], [contour.y]]
     126
     127        #close the contour
     128        if x[-1] != x[0] or y[-1] != y[0]:
     129            x.append(x[0])
     130            y.append(y[0])
     131
     132        setattr(out, 'x', x)
     133        setattr(out, 'y', y)
     134        setattr(out, 'nods', len(x))
     135
     136        return out
     137    #}}}
     138
     139    def checkconsistency(self, *args): #{{{
     140        #recover options
     141        options = pairoptions(*args)
     142
     143        #figure out if two boundaries are identical
     144        for i in range(len(self.boundaries)):
     145            namei = self.boundaries[i].shpfilename
     146            for j in range(1, len(self.boundaries)):
     147                namej = self.boundaries[j].shpfilename
     148                if namei == namej:
     149                    raise RuntimeError('Basin {} has two identical boundaries named {}'.format(self.name, namei))
     150
     151        #go through boundaries and check their consistency
     152        for i in range(len(self.boundaries)):
     153            boundary == self.boundaries[i]
     154            boundary.checkconsistency()
     155    #}}}
     156
     157    def contourplot(self, *args): #{{{
     158        contour = self.contour()
     159        plot(contour.x, contour.y, 'r*-')
     160    #}}}
     161
     162    def shapefilecrop(self, *args): #{{{
     163        #recover options
     164        options = pairoptions(*args)
     165        threshold = options.getfieldvalue('threshold', .65) #.65 degrees lat, long
     166        inshapefile = options.getfieldvalue('shapefile')
     167        outshapefile = options.getfieldvalue('outshapefile', '')
     168
     169        if options.exist('epsgshapefile') and options.exist('projshapefile'):
     170            raise RuntimeError('Basin shapefilecrop error messgae: cannot specify epsgshapefile and projshapefile at the same time')
     171
     172        if options.exist('epsgshapefile'):
     173            projshapefile = epsg2proj(options.getfieldvalue('epsgshapefile'))
     174        elif options.exist('projshapefile'):
     175            projshapefile = options.getfieldvalue('projshapefile')
     176        else:
     177            raise RuntimeError('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified')
     178
     179        #create list of contours that have critical length > threshold (in lat, long)
     180        contours = shpread(inshapefile)
     181        llist = []
     182        for i in range(len(contours)):
     183            contour = contours[i]
     184            carea = polyarea(contour.x, contour.y)
     185            clength = math.sqrt(carea)
     186            if clength < threshold:
     187                llist = np.concatenate(llist, i)
     188        setattr(contours, llist, [])
     189
     190        #project onto reference frame
     191        for i in range(len(contours)):
     192            h = contours[i]
     193            h.x, h.y = gdaltransform(h.x, h.y, projshapefile, self.proj)
     194            contours[i].x = h.x
     195            contours[i].y = h.y
     196
     197        #only keep the contours that are inside the basin (take centroids)
     198        ba = self.contour()
     199        flags = np.zeros(len(contours))
     200        for i in range(len(contours))
     201            h = contours[i]
     202            isin = inpolygon(h.x, h.y, ba.x, ba.y)
     203            if len(np.where(isin == 0, isin, isin)):
     204                flags[i] = 1
     205        pos = flags.nonzero()
     206        contours[pos] = []
     207
     208        #Two options
     209        if outputshapefile == '':
     210            output = contours
     211        else:
     212            shpwrite(contours, outputshapefile)
     213    #}}}
     214#}}}

  • ../trunk-jpl/src/m/classes/boundary.m

     
    1919        end
    2020        methods
    2121                function self = boundary(varargin) % {{{
    22                         switch nargin
    23                                 case 0
    24                                         self=setdefaultparameters(self);
    25                                 otherwise
    26                                         self=setdefaultparameters(self);
    27                                         options=pairoptions(varargin{:});
    28                        
    29                                         %recover field values:
    30                                         self.shppath=getfieldvalue(options,'shppath','');
    31                                         self.shpfilename=getfieldvalue(options,'shpfilename','');
    32                                         self.orientation=getfieldvalue(options,'orientation','normal');
     22                        self=setdefaultparameters(self);
    3323
    34                                         %figure out projection string:
     24                        if nargin > 0,
     25                                options=pairoptions(varargin{:});
     26
     27                                %recover field values:
     28                                self.shppath=getfieldvalue(options,'shppath','');
     29                                self.shpfilename=getfieldvalue(options,'shpfilename','');
     30                                self.orientation=getfieldvalue(options,'orientation','normal');
     31
     32                                %figure out projection string:
     33                                if exist(options,'epsg'),
     34                                        if exist(options,'proj'),
     35                                                error(['Error in constructor for boundary ' self.shppath '. Cannot supply epsg and proj at the same time!']);
     36                                        end
     37                                        epsg=getfieldvalue(options,'epsg');
     38                                        proj=epsg2proj(epsg); % convert to PROJ.4 format
     39                                elseif exist(options,'proj'),
    3540                                        if exist(options,'epsg'),
    36                                                 if exist(options,'proj'),
    37                                                         error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
    38                                                 end
    39                                                 epsg=getfieldvalue(options,'epsg');
    40                                                 %convert into proj4:
    41                                                 proj=epsg2proj(epsg);
    42                                         elseif exist(options,'proj'),
    43                                                 if exist(options,'epsg'),
    44                                                         error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
    45                                                 end
    46                                                 proj=getfieldvalue(options,'proj');
    47                                         else
    48                                                 proj= '+proj=longlat +datum=WGS84 +no_defs';
     41                                                error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']);
    4942                                        end
    50                                         self.proj=proj;
     43                                        proj=getfieldvalue(options,'proj');
     44                                else
     45                                        proj= '+proj=longlat +datum=WGS84 +no_defs';
     46                                end
     47                                self.proj=proj;
    5148                        end
    5249                end % }}}
    5350                function self = setdefaultparameters(self) % {{{
    … …  
    5451                self.shppath='';
    5552                self.shpfilename='';
    5653                self.orientation='normal';
    57                 self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
     54                self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326
    5855                end % }}}
    5956                function disp(self) % {{{
    6057                        disp(sprintf('   boundary parameters:'));
    … …  
    6966                        output=self.shpfilename;
    7067                end % }}}
    7168                function output=edges(self) % {{{
    72                
    73                 %read domain:
    74                 [path,name,ext]=fileparts(self.shpfilename);
    75                 if ~strcmpi(ext,'shp'),
    76                         ext='shp';
    77                 end
    78                 output=shpread([self.shppath '/' name '.' ext]);
    7969
    80                 %do we reverse?
    81                 if strcmpi(self.orientation,'reverse'),
    82                         output.x=flipud(output.x);
    83                         output.y=flipud(output.y);
    84                 end
     70                        %read domain:
     71                        [path,name,ext]=fileparts(self.shpfilename);
     72                        if ~strcmpi(ext,'shp'),
     73                                ext='shp';
     74                        end
     75                        output=shpread([self.shppath '/' name '.' ext]);
     76
     77                        %do we reverse?
     78                        if strcmpi(self.orientation,'reverse'),
     79                                output.x=flipud(output.x);
     80                                output.y=flipud(output.y);
     81                        end
    8582                end % }}}
    8683                function plot(self,varargin) % {{{
    8784                        %recover options
    … …  
    10097                        fontsize=getfieldvalue(options,'fontsize',10);
    10198                        label=getfieldvalue(options,'label','on');
    10299
    103                         if (exist(options,'epsg') &  exist(options,'proj')),
    104                                 error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
    105                         end
    106                         if exist(options,'epsg'),
     100                        if (exist(options,'epsg'),
     101                                if exist(options,'proj')),
     102                                        error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
    107103                                proj=epsg2proj(getfieldvalue(options,'epsg'));
    108104                        elseif exist(options,'proj'),
    109105                                proj=getfieldvalue(options,'proj');
    … …  
    118114                        end
    119115                        domain=shpread([self.shppath '/' name '.' ext]);
    120116
    121                         %convert boundary to another reference frame:  {{{
     117                        %convert boundary to another reference frame: {{{
    122118                        for i=1:length(domain),
    123119                                try,
    124120                                        [x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj);
    … …  
    168164                                dmax=max(distance);
    169165                                isidentical=find(distance<dmax*tolerance);
    170166                                if ~isempty(isidentical),
    171                                         error(['boundary ' shpfilename ' has two  vertices extremely close to one another']);
     167                                        error(['boundary ' shpfilename ' has two vertices extremely close to one another']);
    172168                                end
    173169                        end
    174170                end % }}}
    … …  
    189185                        fontsize=getfieldvalue(options,'fontsize',10);
    190186                        label=getfieldvalue(options,'label','on');
    191187
     188                        if (exist(options,'epsg'),
     189                                if exist(options,'proj')),
     190                                        error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
     191                                proj=epsg2proj(getfieldvalue(options,'epsg'));
     192                        elseif exist(options,'proj'),
     193                                proj=getfieldvalue(options,'proj');
     194                        else
     195                                proj=epsg2proj(4326);
     196                        end
     197
    192198                        %read domain:
    193199                        [path,name,ext]=fileparts(self.shpfilename);
    194200                        if ~strcmpi(ext,'shp'),

  • ../trunk-jpl/src/m/plot/plotmodel.py

     
    1 import numpy as np
    2 from plotoptions import plotoptions
    3 from plot_manager import plot_manager
    41from math import ceil, sqrt
    52
    63try:
    … …  
    85    from mpl_toolkits.axes_grid1 import ImageGrid
    96except ImportError:
    107    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     8import numpy as np
    119
     10from plot_manager import plot_manager
     11from plotoptions import plotoptions
    1212
     13
    1314def plotmodel(md, *args):
    14     ''' at command prompt, type 'plotdoc()' for additional documentation
    1515    '''
     16    PLOTMODEL - At command prompt, type 'plotdoc()' for additional
     17    documentation.
    1618
     19    TODO:
     20        - Fix 'plotdoc()', as it is not currently working.
     21    '''
     22
    1723    #First process options
    1824    options = plotoptions(*args)
    19     #get number of subplots
    20     subplotwidth = ceil(sqrt(options.numberofplots))
     25
    2126    #Get figure number and number of plots
    2227    figurenumber = options.figurenumber
    2328    numberofplots = options.numberofplots
    2429
    25     #get hold
    26     hold = options.list[0].getfieldvalue('hold', False)
     30    #get number of subplots
     31    subplotwidth = ceil(sqrt(numberofplots))
    2732
     33    # TODO: Check that commenting this out is correct; we do not need a hold
     34    # under matplotlib, right?
     35    #
     36    # #get hold
     37    # hold = options.list[0].getfieldvalue('hold', False)
     38
    2839    #if nrows and ncols specified, then bypass
    2940    if options.list[0].exist('nrows'):
    3041        nrows = options.list[0].getfieldvalue('nrows')
    … …  
    4354    nrows = int(nrows)
    4455
    4556    #check that nrows and ncols were given at the same time!
    46     if not nr == nc:
    47         raise Exception('error: nrows and ncols need to be specified together, or not at all')
     57    if nr != nc:
     58        raise Exception('plotmodel error message: nrows and ncols need to be specified together, or not at all')
    4859
    49     #Go through plots
     60    # Go through plots
     61    #
     62    # NOTE: The following is where Python + matplolib differs substantially in
     63    #       implementation and inteface from MATLAB.
     64    #
     65    # Sources:
     66    # - https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html
     67    #
    5068    if numberofplots:
    51         #if plt.fignum_exists(figurenumber):
    52         #       plt.cla()
    53 
    5469        #if figsize specified
    5570        if options.list[0].exist('figsize'):
    5671            figsize = options.list[0].getfieldvalue('figsize')
    … …  
    6277        backgroundcolor = options.list[0].getfieldvalue('backgroundcolor', (0.7, 0.7, 0.7))
    6378        fig.set_facecolor(backgroundcolor)
    6479
    65         translator = {'on': 'each',
    66                       'off': 'None',
    67                       'one': 'single'}
    68     # options needed to define plot grid
     80        # options needed to define plot grid
    6981        plotnum = options.numberofplots
    7082        if plotnum == 1:
    7183            plotnum = None
    72         direction = options.list[0].getfieldvalue('direction', 'row')  # row, column
    73         axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25)
    74         add_all = options.list[0].getfieldvalue('add_all', True)  # True, False
    75         share_all = options.list[0].getfieldvalue('share_all', True)  # True, False
    76         label_mode = options.list[0].getfieldvalue('label_mode', 'L')  # 1, L, all
     84
     85        # NOTE: The inline comments for each of the following parameters are
     86        #       taken from https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html
     87        #
     88        direction = options.list[0].getfieldvalue('direction', 'row') # {"row", "column"}, default: "row"
     89        axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25) # float or (float, float), default : 0.02; Padding or (horizonal padding, vertical padding) between axes, in inches
     90        add_all = options.list[0].getfieldvalue('add_all', True) # bool, default: True
     91        share_all = options.list[0].getfieldvalue('share_all', True) # bool, default: False
     92        label_mode = options.list[0].getfieldvalue('label_mode', 'L') # {"L", "1", "all"}, default: "L"; Determines which axes will get tick labels: "L": All axes on the left column get vertical tick labels; all axes on the bottom row get horizontal tick labels;. "1": Only the bottom left axes is labelled. "all": all axes are labelled.
     93
     94        # Translate MATLAB colorbar mode to matplotlib
     95        #
     96        # TODO:
     97        # - Add 'edge' option (research if there is a corresponding option in
     98        #   MATLAB)?
     99        #
    77100        colorbar = options.list[0].getfieldvalue('colorbar', 'on')  # on, off (single)
    78         cbar_mode = translator[colorbar]
    79         cbar_location = options.list[0].getfieldvalue('colorbarpos', 'right')  # right, top
    80         cbar_size = options.list[0].getfieldvalue('colorbarsize', '5%')
    81         cbar_pad = options.list[0].getfieldvalue('colorbarpad', 0.025)  # None or %
    82101
    83         axgrid = ImageGrid(fig, 111,
    84                            nrows_ncols=(nrows, ncols),
    85                            direction=direction,
    86                            axes_pad=axes_pad,
    87                            add_all=add_all,
    88                            share_all=share_all,
    89                            label_mode=label_mode,
    90                            cbar_mode=cbar_mode,
    91                            cbar_location=cbar_location,
    92                            cbar_size=cbar_size,
    93                            cbar_pad=cbar_pad)
     102        if colorbar == 'on':
     103            colorbar = 'each'
     104        elif colorbar == 'one':
     105            colorbar = 'single'
     106        elif colorbar == 'off':
     107            colorbar = 'None'
     108        else:
     109            raise RuntimeError('plotmodel error: colorbar mode \'{}\' is not a valid option'.format(colorbar))
    94110
     111        cbar_mode = colorbar # {"each", "single", "edge", None }, default: None
     112        cbar_location = options.list[0].getfieldvalue('colorbarpos', 'right') # {"left", "right", "bottom", "top"}, default: "right"
     113        cbar_pad = options.list[0].getfieldvalue('colorbarpad', 0.025) # float, default: None
     114        cbar_size = options.list[0].getfieldvalue('colorbarsize', '5%') # size specification (see Size.from_any), default: "5%"
     115
     116        # NOTE: Second parameter is:
     117        #
     118        #   rect(float, float, float, float) or int
     119        #
     120        # The axes position, as a (left, bottom, width, height) tuple or as a
     121        # three-digit subplot position code (e.g., "121").
     122        #
     123        axgrid = ImageGrid(
     124            fig,
     125            111,
     126            nrows_ncols=(nrows, ncols),
     127            direction=direction,
     128            axes_pad=axes_pad,
     129            add_all=add_all,
     130            share_all=share_all,
     131            label_mode=label_mode,
     132            cbar_mode=cbar_mode,
     133            cbar_location=cbar_location,
     134            cbar_size=cbar_size,
     135            cbar_pad=cbar_pad
     136        )
     137
    95138        if cbar_mode == 'None':
    96139            for ax in axgrid.cbar_axes:
    97140                fig._axstack.remove(ax)

  • ../trunk-jpl/src/m/plot/plot_unit.py

     
    11from cmaptools import getcolormap, truncate_colormap
    2 from plot_quiver import plot_quiver
    3 import numpy as np
    42try:
    53    import matplotlib as mpl
    64    import matplotlib.pyplot as plt
    … …  
    97    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    108except ImportError:
    119    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     10import numpy as np
    1211
     12from plot_quiver import plot_quiver
    1313
     14
    1415def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, axgrid, gridindex):
    15     """
     16    '''
    1617    PLOT_UNIT - unit plot, display data
    1718
    1819    Usage:
    19     plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options)
     20        plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options)
    2021
    2122    See also: PLOTMODEL, PLOT_MANAGER
    22     """
     23    '''
    2324    #if we are plotting 3d replace the current axis
    2425    if not is2d:
    2526        axgrid[gridindex].axis('off')
    … …  
    227228        else:
    228229            raise ValueError('plot_unit error: 3D node plot not supported yet')
    229230        return
    230 
    231231    # }}}
    232232    # {{{ plotting P1 Patch (TODO)
    233 
    234233    elif datatype == 4:
    235234        print('plot_unit message: P1 patch plot not implemented yet')
    236235        return
    237 
    238236    # }}}
    239237    # {{{ plotting P0 Patch (TODO)
    240 
    241238    elif datatype == 5:
    242239        print('plot_unit message: P0 patch plot not implemented yet')
    243240        return
    244 
    245241    # }}}
    246242    else:
    247243        raise ValueError('datatype = %d not supported' % datatype)

  • ../trunk-jpl/src/m/plot/applyoptions.py

     
    1 import numpy as np
    21from cmaptools import getcolormap
    3 from plot_contour import plot_contour
    4 from plot_streamlines import plot_streamlines
    5 from expdisp import expdisp
    6 
    72try:
    8     from matplotlib.ticker import MaxNLocator
    93    import matplotlib as mpl
    104    import matplotlib.pyplot as plt
     5    from matplotlib.ticker import MaxNLocator
    116except ImportError:
    127    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     8import numpy as np
    139
     10from expdisp import expdisp
     11from plot_contour import plot_contour
     12from plot_streamlines import plot_streamlines
    1413
     14
    1515def applyoptions(md, data, options, fig, axgrid, gridindex):
    1616    '''
    1717    APPLYOPTIONS - apply options to current plot
    … …  
    1919    'plotobj' is the object returned by the specific plot call used to
    2020    render the data.  This object is used for adding a colorbar.
    2121
    22         Usage:
    23             applyoptions(md, data, options)
     22    Usage:
     23        applyoptions(md, data, options)
    2424
    25         See also: PLOTMODEL, PARSE_OPTIONS
     25    See also: PLOTMODEL, PARSE_OPTIONS
    2626    '''
    2727
    2828    # get handle to current figure and axes instance
    … …  
    3333    fontsize = options.getfieldvalue('fontsize', 8)
    3434    fontweight = options.getfieldvalue('fontweight', 'normal')
    3535    fontfamily = options.getfieldvalue('fontfamily', 'sans - serif')
    36     font = {'fontsize': fontsize,
    37             'fontweight': fontweight,
    38             'family': fontfamily}
     36    font = {
     37        'fontsize': fontsize,
     38        'fontweight': fontweight,
     39        'family': fontfamily
     40    }
    3941    # }}}
    4042    # {{{ title
    4143    if options.exist('title'):
Note: See TracBrowser for help on using the repository browser.