Changeset 25065


Ignore:
Timestamp:
06/18/20 16:01:39 (5 years ago)
Author:
jdquinn
Message:

CHG: Saving progress on SLR/Solid Earth translation

Location:
issm/trunk-jpl/src/m
Files:
7 added
21 edited

Legend:

Unmodified
Added
Removed

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

    r24779 r25065  
    1919        methods
    2020                function self = basin(varargin) % {{{
    21                         switch nargin
    22                                 case 0
    23                                         self=setdefaultparameters(self);
    24                                 otherwise
    25 
    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','');
    33 
    34                                         %figure out projection string:
     21                        self=setdefaultparameters(self);
     22
     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','');
     30
     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 % }}}
    … …  
    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
    … …  
    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 % }}}
    … …  
    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
    … …  
    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=[];
    … …  
    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

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

    r24779 r25065  
    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');
    33 
    34                                         %figure out projection string:
     22                        self=setdefaultparameters(self);
     23
     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';
    49                                         end
    50                                         self.proj=proj;
     41                                                error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']);
     42                                        end
     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 % }}}
    … …  
    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) % {{{
    … …  
    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]);
    79 
    80                 %do we reverse?
    81                 if strcmpi(self.orientation,'reverse'),
    82                         output.x=flipud(output.x);
    83                         output.y=flipud(output.y);
    84                 end
     69
     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) % {{{
    … …  
    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'),
    … …  
    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,
    … …  
    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
    … …  
    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);

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

    r25012 r25065  
    11from collections import OrderedDict
     2
    23import pairoptions
    34
    … …  
    78    PLOTOPTIONS class definition
    89
    9         Usage:
    10             plotoptions = plotoptions(*arg)
     10    Usage:
     11        plotoptions = plotoptions(*arg)
    1112    '''
    1213
    … …  
    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)
    … …  
    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])

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

    r24886 r25065  
    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) % {{{

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

    r25035 r25065  
    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
    … …  
    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
    … …  
    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:
    … …  
    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
    … …  
    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
    … …  
    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
    … …  
    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))

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

    r25035 r25065  
    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

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

    r25035 r25065  
    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
    … …  
    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
    … …  
    5357    os.remove(filename_in)
    5458    os.remove(filename_out)
     59
     60    return [xout, yout]
    5561#}}}

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

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

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

    r25035 r25065  
    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

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

    r24942 r25065  
    1 import os
    21import subprocess
    32
    … …  
    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
    … …  
    3433
    3534    #process options
    36     options = pairoptions(*varargin)
     35    options = pairoptions(*args)
    3736    #options = deleteduplicates(options, 1)
    3837

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

    r24745 r25065  
    307307        end
    308308end
    309 
    310309
    311310%text (default value is empty, not NaN...)
    … …  
    352351end
    353352
    354 %streamliness
     353%streamlines
    355354if exist(options,'streamlines'),
    356355        plot_streamlines(md,options);

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

    r24565 r25065  
     1from cmaptools import getcolormap
     2try:
     3    import matplotlib as mpl
     4    import matplotlib.pyplot as plt
     5    from matplotlib.ticker import MaxNLocator
     6except ImportError:
     7    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
    18import numpy as np
    2 from cmaptools import getcolormap
     9
     10from expdisp import expdisp
    311from plot_contour import plot_contour
    412from plot_streamlines import plot_streamlines
    5 from expdisp import expdisp
    6 
    7 try:
    8     from matplotlib.ticker import MaxNLocator
    9     import matplotlib as mpl
    10     import matplotlib.pyplot as plt
    11 except ImportError:
    12     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
    1313
    1414
    … …  
    2020    render the data.  This object is used for adding a colorbar.
    2121
    22         Usage:
    23             applyoptions(md, data, options)
    24 
    25         See also: PLOTMODEL, PARSE_OPTIONS
     22    Usage:
     23        applyoptions(md, data, options)
     24
     25    See also: PLOTMODEL, PARSE_OPTIONS
    2626    '''
    2727
    … …  
    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

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

    r25038 r25065  
    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);
    … …  
    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);
    … …  
    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',...
    … …  
    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 ';']);
    … …  
    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);
    … …  
    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);
    … …  
    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);

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

    r24269 r25065  
    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
    … …  
    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
    … …  
    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.
     29
     30    'gridindex' is passed through to specialized plot* functions and
     31    applyoptions.
    3132
    3233    Usage:
    … …  
    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
    … …  
    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':

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

    r24256 r25065  
    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
    … …  
    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
     18
    1719
    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
    … …  
    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

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

    r24569 r25065  
    11from cmaptools import getcolormap, truncate_colormap
    2 from plot_quiver import plot_quiver
    3 import numpy as np
    42try:
    53    import matplotlib as mpl
    … …  
    108except ImportError:
    119    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     10import numpy as np
     11
     12from plot_quiver import plot_quiver
    1213
    1314
    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:
    … …  
    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:

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

    r25038 r25065  
    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.
    … …  
    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
    … …  
    125125        end
    126126else
    127         error('plotmodel error message: no output data found. ');
     127        error('plotmodel error message: no output data found.');
    128128end

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

    r24269 r25065  
    1 import numpy as np
    2 from plotoptions import plotoptions
    3 from plot_manager import plot_manager
    41from math import ceil, sqrt
    52
    … …  
    96except ImportError:
    107    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     8import numpy as np
     9
     10from plot_manager import plot_manager
     11from plotoptions import plotoptions
    1112
    1213
    1314def plotmodel(md, *args):
    14     ''' at command prompt, type 'plotdoc()' for additional documentation
     15    '''
     16    PLOTMODEL - At command prompt, type 'plotdoc()' for additional
     17    documentation.
     18
     19    TODO:
     20        - Fix 'plotdoc()', as it is not currently working.
    1521    '''
    1622
    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))
     32
     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)
    2738
    2839    #if nrows and ncols specified, then bypass
    … …  
    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'):
    … …  
    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))
     110
     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        )
    94137
    95138        if cbar_mode == 'None':

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

    r24894 r25065  
    1 import numpy as np
    21from collections import OrderedDict
    32from copy import deepcopy
     3
     4import numpy as np
    45
    56
    … …  
    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']
    15 
    16 Note that 'x' returns the array and x.__dict__ will only return
    17     attributes other than the array
    18 
    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
    23 
    24 Example uses:
    25 
    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
    38 
    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')
    44 
    45 Credit: https://github.com/Vectorized/Python-Attribute-List
    46 '''
     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']
     16
     17    Note that 'x' returns the array and x.__dict__ will only return attributes
     18    other than the array
     19
     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
     24
     25    Examples:
     26
     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
     39
     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')
     45
     46    Sources:
     47    -https://github.com/Vectorized/Python-Attribute-List
     48    '''
    4749
    4850    def __new__(self, *args, **kwargs):
    … …  
    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.
    68 
    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
     67    A form of dictionary-like structure that maintains the ordering in which
     68    its fields/attributes and their corresponding values were added.
     69
     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
    77 
    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
    82 
    83 Supports: len(x), str(x), for-loop iteration.
    84 Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
    85 
    86 Usage:
    87     x = OrderedStruct()
    88     x.y = 5
    89     x.z = 6
    90     OR
    91     x = OrderedStruct('y', 5, 'z', 6)
    92 
    93     # note below that the output fields as iterables are always
    94     #    in the same order as the inputs
     74    Example:
     75        OrderedDict:  # a bit clumsy to use and look at
     76            x['y'] = 5
     77
     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
     82
     83    Supports: len(x), str(x), for-loop iteration.
     84    Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
     85
     86    Usage:
     87        x = OrderedStruct()
     88        x.y = 5
     89        x.z = 6
     90        OR
     91        x = OrderedStruct('y', 5, 'z', 6)
     92
     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']
    … …  
    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'''
    121 
    122     # keys and values
     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        '''
     123
     124        # keys and values
    123125        self._k = []
    124126        self._v = []
    … …  
    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()):
    … …  
    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()):
    … …  
    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
    … …  
    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]:
    … …  
    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
    … …  
    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
    … …  
    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
    … …  
    297312def fullfile(*args):
    298313    '''
    299     use:
    300 
    301     fullfile(path, path, ... , file_name + ext)
     314    usage:
     315        fullfile(path, path, ... , file_name + ext)
     316
    302317    returns: "path/path/.../file_name.ext"
    303318
    … …  
    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:])):
    … …  
    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:
    … …  
    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
    339 
    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
    342 
    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 '''
     356    as_numpy_ndarray will return the result as a numpy.ndarray and is False by
     357    default
     358
     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
     361
     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)

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

    r22921 r25065  
    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
    … …  
    2727end
    2828
    29 %initialize number of profile
    30 count=0;
    31 
    3229%read shapefile
    3330shp=shaperead(filename);
    … …  
    3633fields=fieldnames(shp);
    3734for i=1:length(shp),
    38         if strcmpi(shp(i).Geometry,'Polygon'),
     35        if strcmpi(shp(i).Geometry,'Point'),
     36                x=shp(i).X'; y=shp(i).Y';
     37                ids=find(isnan(x));
     38                x(ids)=[]; y(ids)=[];
     39
     40                Struct(end+1).x=x;
     41                Struct(end).y=y;
     42                Struct(end).density=1;
     43                if isfield(shp,'id'),
     44                        Struct(end).name=num2str(shp(i).id);
     45                else
     46                        Struct(end).name='';
     47                end
     48                for j=1:length(fields),
     49                        field=fields{j};
     50                        if ~(strcmpi(field,'X') | strcmpi(field,'Y') | strcmpi(field,'id')),
     51                                Struct(end).(field)=shp(i).(field);
     52                        end
     53                end
     54        elseif strcmpi(shp(i).Geometry,'Line'),
    3955                x=shp(i).X'; y=shp(i).Y';
    4056                ids=find(isnan(x));
    … …  
    5773                        end
    5874                end
    59         end
    60        
    61         if strcmpi(shp(i).Geometry,'Line'),
     75        elseif strcmpi(shp(i).Geometry,'Polygon'),
    6276                x=shp(i).X'; y=shp(i).Y';
    6377                ids=find(isnan(x));
    … …  
    6983                Struct(end).density=1;
    7084                Struct(end).closed=1;
    71                 if isfield(shp,'id'),
    72                         Struct(end).name=num2str(shp(i).id);
    73                 else
    74                         Struct(end).name='';
    75                 end
    76                 for j=1:length(fields),
    77                         field=fields{j};
    78                         if ~(strcmpi(field,'X') | strcmpi(field,'Y') | strcmpi(field,'id')),
    79                                 Struct(end).(field)=shp(i).(field);
    80                         end
    81                 end
    82         end
    83 
    84 
    85         if strcmpi(shp(i).Geometry,'Point'),
    86                 x=shp(i).X'; y=shp(i).Y';
    87                 ids=find(isnan(x));
    88                 x(ids)=[]; y(ids)=[];
    89 
    90                 Struct(end+1).x=x;
    91                 Struct(end).y=y;
    92                 Struct(end).density=1;
    9385                if isfield(shp,'id'),
    9486                        Struct(end).name=num2str(shp(i).id);

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

    r22920 r25065  
    1010%   See also SHPREAD
    1111
    12 
    13 %initialize number of profile
    14 count=0;
    15 
    1612contours=struct([]);
    1713for i=1:length(shp),
    … …  
    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;
    … …  
    2725        contours(i).Y=shp(i).y;
    2826end
    29        
     27
    3028shapewrite(contours,filename);
    3129end
Note: See TracChangeset for help on using the changeset viewer.