Index: ../trunk-jpl/src/m/shp/shpwrite.py =================================================================== --- ../trunk-jpl/src/m/shp/shpwrite.py (nonexistent) +++ ../trunk-jpl/src/m/shp/shpwrite.py (revision 25065) @@ -0,0 +1,60 @@ +try: + import shapefile +except ImportError: + print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled") + +def shpwrite(shp, filename): #{{{ + ''' + SHPREAD - write a shape file from a contour structure + + The current implementation of shpwrite depends on PyShp. + + Usage: + shpwrite(shp, filename) + + Example: + shpwrite(shp, 'domainoutline.shp') + + See also SHPREAD + + Sources: + - https://github.com/GeospatialPython/pyshp + + TODO: + - Should we check if there is only one element (see how MATLAB's shaperead + and shapewrite handle single shapes versus multiple shapes)? + ''' + + # Sanity checks + for shape in shp: + print(shape) + + if shp[0].Geometry == 'Point': + shapeType = 1 + elif shp[0].Geometry == 'Line': + shapeType = 3 + elif shp[0].Geometry == 'Polygon': + shapeType = 5 + else: + raise Exception('shpwrite error: geometry \'{}\' is not currently supported'.format(shp[0].Geometry)) + + sf = shapefile.Writer(filename, shapeType=shapeType) + + for i in range(len(shp)): + sf.field('name', 'C') # TODO: Allow shape ids/names to be passed through + if shapeType == 1: # POINT + sf.point(shp[i].x, shp[i].y) + elif shapeType == 3: # POLYLINE + points = [] + for j in range(len(shp[i].x)): + points.append([shp[i].x[j], shp[i].y[j]]) + sf.line(line) + elif shapeType == 5: # POLYGON + points = [] + for j in range(len(shp[i].x)): + points.append([shp[i].x[j], shp[i].y[j]]) + sf.poly(points) + sf.null() + sf.record(str(i)) + sf.close() +#}}} Index: ../trunk-jpl/src/m/shp/shpwrite.m =================================================================== --- ../trunk-jpl/src/m/shp/shpwrite.m (revision 25064) +++ ../trunk-jpl/src/m/shp/shpwrite.m (revision 25065) @@ -9,10 +9,6 @@ % % See also SHPREAD - -%initialize number of profile -count=0; - contours=struct([]); for i=1:length(shp), if strcmpi(shp(i).Geometry,'Point'), @@ -21,11 +17,13 @@ contours(i).Geometry='Polygon'; elseif strcmpi(shp(i).Geometry,'Line'), contours(i).Geometry='Line'; + else, + error(['shpwrite error: geometry ' shp(i).Geometry ' not currently supported']); end contours(i).id=i; contours(i).X=shp(i).x; contours(i).Y=shp(i).y; end - + shapewrite(contours,filename); end Index: ../trunk-jpl/src/m/mesh/meshintersect3d.py =================================================================== --- ../trunk-jpl/src/m/mesh/meshintersect3d.py (revision 25064) +++ ../trunk-jpl/src/m/mesh/meshintersect3d.py (revision 25065) @@ -6,9 +6,11 @@ from pairoptions import pairoptions + def meshintersect3d(x, y, z, xs, ys, zs, *args): #{{{ ''' - 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). + 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). ''' # Process options Index: ../trunk-jpl/src/m/coordsystems/laea.py =================================================================== --- ../trunk-jpl/src/m/coordsystems/laea.py (nonexistent) +++ ../trunk-jpl/src/m/coordsystems/laea.py (revision 25065) @@ -0,0 +1,15 @@ +def laea(lat, long): #{{{ + ''' + LAEA - Lambert Azimuthal Equal Area projection at lat, long projection + center. + + Usage: + string = laea(45, -90) + + Example: + string = laea(45, -90) + return string = '+proj=laea +lat_0=45 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs' + ''' + + return '+proj=laea +lat_0={} +lon_0={} +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'.format(lat, long) +#}}} Index: ../trunk-jpl/src/m/coordsystems/gdaltransform.m =================================================================== --- ../trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25064) +++ ../trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25065) @@ -19,9 +19,9 @@ % Bamber's Greeland projection % +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 % -% To get proj.4 string from EPSG, use gdalsrsinfo. Example: +% To get proj.4 string from EPSG, use gdalsrsinfo. Example: % -% gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" +% gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" %give ourselves unique file names filename_in = tempname(); Index: ../trunk-jpl/src/m/geometry/AboveGround.py =================================================================== --- ../trunk-jpl/src/m/geometry/AboveGround.py (nonexistent) +++ ../trunk-jpl/src/m/geometry/AboveGround.py (revision 25065) @@ -0,0 +1,8 @@ +import numpy as np + +def function AboveGround(lat, long, r, height): #{{{ + r = r + height + x = r * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(long)) + y = r * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(long)) + z = r * np.sin(np.deg2rad(lat)) +#}}} Index: ../trunk-jpl/src/m/classes/plotoptions.py =================================================================== --- ../trunk-jpl/src/m/classes/plotoptions.py (revision 25064) +++ ../trunk-jpl/src/m/classes/plotoptions.py (revision 25065) @@ -1,4 +1,5 @@ from collections import OrderedDict + import pairoptions @@ -6,8 +7,8 @@ ''' PLOTOPTIONS class definition - Usage: - plotoptions = plotoptions(*arg) + Usage: + plotoptions = plotoptions(*arg) ''' def __init__(self, *arg): # {{{ @@ -35,7 +36,7 @@ def buildlist(self, *arg): #{{{ #check length of input if len(arg) % 2: - raise TypeError('Invalid parameter / value pair arguments') + raise TypeError('Invalid parameter/value pair arguments') #go through args and build list (like pairoptions) rawoptions = pairoptions.pairoptions(*arg) @@ -102,7 +103,7 @@ if len(nums) != 2: continue if False in [x.isdigit() for x in nums]: - raise ValueError('error: in option i - j both i and j must be integers') + raise ValueError('error: in option i-j both i and j must be integers') for j in range(int(nums[0]) - 1, int(nums[1])): self.list[j].addfield(field, rawlist[i][1]) Index: ../trunk-jpl/src/m/classes/basin.m =================================================================== --- ../trunk-jpl/src/m/classes/basin.m (revision 25064) +++ ../trunk-jpl/src/m/classes/basin.m (revision 25065) @@ -18,42 +18,39 @@ end methods function self = basin(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - otherwise + self=setdefaultparameters(self); - self=setdefaultparameters(self); - options=pairoptions(varargin{:}); - - %recover field values: - self.boundaries=getfieldvalue(options,'boundaries',{}); - self.name=getfieldvalue(options,'name',''); - self.continent=getfieldvalue(options,'continent',''); + if nargin>0, + options=pairoptions(varargin{:}); + + %recover field values: + self.boundaries=getfieldvalue(options,'boundaries',{}); + self.name=getfieldvalue(options,'name',''); + self.continent=getfieldvalue(options,'continent',''); - %figure out projection string: + %figure out projection string: + if exist(options,'epsg'), + if exist(options,'proj'), + error(['Error in constructor for basin ' self.name '. Cannot supply epsg and proj at the same time!']); + end + epsg=getfieldvalue(options,'epsg'); + proj=epsg2proj(epsg); %convert to PROJ.4 format + elseif exist(options,'proj'), if exist(options,'epsg'), - if exist(options,'proj'), - error(['Error in constructor for boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']); - end - epsg=getfieldvalue(options,'epsg'); - %convert into proj4: - proj=epsg2proj(epsg); - elseif exist(options,'proj'), - if exist(options,'epsg'), - error(['Error in constructor for boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']); - end - proj=getfieldvalue(options,'proj'); - else - proj= '+proj=longlat +datum=WGS84 +no_defs'; + error(['Error in constructor for basin ' self.name '. Cannot supply proj and epsg at the same time!']); end - self.proj=proj; + proj=getfieldvalue(options,'proj'); + else + proj='+proj=longlat +datum=WGS84 +no_defs'; + end + + self.proj=proj; end end % }}} function self = setdefaultparameters(self) % {{{ self.name=''; self.continent=''; - self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326 + self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326 self.boundaries={}; end % }}} @@ -92,11 +89,11 @@ options=pairoptions(varargin{:}); extension=getfieldvalue(options,'extension',1); - [path,nme,ext]=fileparts(self.name); + [path,name,ext]=fileparts(self.name); if extension, - output=[nme ext]; + output=[name ext]; else - output=nme; + output=name; end end % }}} function plot(self,varargin) % {{{ @@ -173,7 +170,7 @@ inshapefile=getfieldvalue(options,'shapefile'); outputshapefile=getfieldvalue(options,'outputshapefile',''); - if (exist(options,'epsgshapefile') & exist(options,'projshapefile')), + if (exist(options,'epsgshapefile') & exist(options,'projshapefile')), error('Basin shapefilecrop error message: cannot specify epsgshapefile and projshapefile at the same time'); end if exist(options,'epsgshapefile'), @@ -185,7 +182,7 @@ end - %create list of contours that have critical length > threshold: (in lat,long) + %create list of contours that have critical length > threshold (in lat,long): contours=shpread(inshapefile); llist=[]; for i=1:length(contours), @@ -211,8 +208,8 @@ flags=zeros(length(contours),1); for i=1:length(contours), h=contours(i); - in=inpolygon(h.x,h.y,ba.x,ba.y); - if ~isempty(find(in==0)), + isin=inpolygon(h.x,h.y,ba.x,ba.y); + if ~isempty(find(isin==0)), flags(i)=1; end end Index: ../trunk-jpl/src/m/qmu/helpers.py =================================================================== --- ../trunk-jpl/src/m/qmu/helpers.py (revision 25064) +++ ../trunk-jpl/src/m/qmu/helpers.py (revision 25065) @@ -1,8 +1,9 @@ -import numpy as np from collections import OrderedDict from copy import deepcopy +import numpy as np + class struct(object): '''An empty struct that can be assigned arbitrary attributes''' pass @@ -10,40 +11,41 @@ class Lstruct(list): ''' -An empty struct that can be assigned arbitrary attributes - but can also be accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd'] + An empty struct that can be assigned arbitrary attributes but can also be + accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd'] -Note that 'x' returns the array and x.__dict__ will only return - attributes other than the array + Note that 'x' returns the array and x.__dict__ will only return attributes + other than the array -List-based and struct-based behaviors work normally, however they are referenced - as if the other does not exist; len(x) corresponds only to - the list component of x, len(x.a) corresponds to x.a, x.__dict__ - corresponds only to the non-x-list attributes + List-based and struct-based behaviors work normally, however they are + referenced as if the other does not exist; len(x) corresponds only to the + list component of x, len(x.a) corresponds to x.a, x.__dict__ corresponds + only to the non-x-list attributes -Example uses: + Examples: -x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4] - x.a = 'hello' - len(x) -> 4 - x.append(5) - len(x) -> 5 - x[2] -> 3 - x.a -> 'hello' - print x -> [1, 2, 3, 4, 5] - x.__dict__ -> {'a':'hello'} - x.b = [6, 7, 8, 9] - x.b[-1] -> 9 - len(x.b) -> 4 + x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4] + x.a = 'hello' + len(x) -> 4 + x.append(5) + len(x) -> 5 + x[2] -> 3 + x.a -> 'hello' + print x -> [1, 2, 3, 4, 5] + x.__dict__ -> {'a':'hello'} + x.b = [6, 7, 8, 9] + x.b[-1] -> 9 + len(x.b) -> 4 -Other valid constructors: - x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3] - x = Lstruct(1, 2, 3)(a = 'hello') - x = Lstruct([1, 2, 3], x = 'hello') - x = Lstruct((1, 2, 3), a = 'hello') + Other valid constructors: + x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3] + x = Lstruct(1, 2, 3)(a = 'hello') + x = Lstruct([1, 2, 3], x = 'hello') + x = Lstruct((1, 2, 3), a = 'hello') -Credit: https://github.com/Vectorized/Python-Attribute-List -''' + Sources: + -https://github.com/Vectorized/Python-Attribute-List + ''' def __new__(self, *args, **kwargs): return super(Lstruct, self).__new__(self, args, kwargs) @@ -62,36 +64,34 @@ class OrderedStruct(object): ''' -A form of dictionary-like structure that maintains the - ordering in which its fields/attributes and their - corresponding values were added. + A form of dictionary-like structure that maintains the ordering in which + its fields/attributes and their corresponding values were added. -OrderedDict is a similar device, however this class - can be used as an "ordered struct/class" giving - it much more flexibility in practice. It is + OrderedDict is a similar device, however this class can be used as an + "ordered struct/class" giving it much more flexibility in practice. It is also easier to work with fixed valued keys in-code. -Eg: -OrderedDict: # a bit clumsy to use and look at - x['y'] = 5 + Example: + OrderedDict: # a bit clumsy to use and look at + x['y'] = 5 -OrderedStruct: # nicer to look at, and works the same way - x.y = 5 - OR - x['y'] = 5 # supports OrderedDict-style usage + OrderedStruct: # nicer to look at, and works the same way + x.y = 5 + OR + x['y'] = 5 # supports OrderedDict-style usage -Supports: len(x), str(x), for-loop iteration. -Has methods: x.keys(), x.values(), x.items(), x.iterkeys() + Supports: len(x), str(x), for-loop iteration. + Has methods: x.keys(), x.values(), x.items(), x.iterkeys() -Usage: - x = OrderedStruct() - x.y = 5 - x.z = 6 - OR - x = OrderedStruct('y', 5, 'z', 6) + Usage: + x = OrderedStruct() + x.y = 5 + x.z = 6 + OR + x = OrderedStruct('y', 5, 'z', 6) - # note below that the output fields as iterables are always - # in the same order as the inputs + note below that the output fields as iterables are always in the same + order as the inputs x.keys() -> ['y', 'z'] x.values() -> [5, 6] @@ -110,16 +110,18 @@ ('x', 5) ('y', 6) - Note: to access internal fields use dir(x) - (input fields will be included, but - are not technically internals) + Note: to access internal fields use dir(x) (input fields will be included, + but are not technically internals) ''' def __init__(self, *args): - '''Provided either nothing or a series of strings, construct a structure that will, -when accessed as a list, return its fields in the same order in which they were provided''' + ''' + Provided either nothing or a series of strings, construct a structure + that will, when accessed as a list, return its fields in the same order + in which they were provided + ''' - # keys and values + # keys and values self._k = [] self._v = [] @@ -184,9 +186,11 @@ yield(a, b) def __copy__(self): - # shallow copy, hard copies of trivial attributes, - # references to structures like lists/OrderedDicts - # unless redefined as an entirely different structure + ''' + shallow copy, hard copies of trivial attributes, + references to structures like lists/OrderedDicts + unless redefined as an entirely different structure + ''' newInstance = type(self)() for k, v in list(self.items()): exec(('newInstance.%s = v') % (k)) @@ -193,11 +197,13 @@ return newInstance def __deepcopy__(self, memo=None): - # hard copy of all attributes - # same thing but call deepcopy recursively - # technically not how it should be done, - # (see https://docs.python.org/2/library/copy.html #copy.deepcopy ) - # but will generally work in this case + ''' + hard copy of all attributes + same thing but call deepcopy recursively + technically not how it should be done, + (see https://docs.python.org/2/library/copy.html #copy.deepcopy ) + but will generally work in this case + ''' newInstance = type(self)() for k, v in list(self.items()): exec(('newInstance.%s = deepcopy(v)') % (k)) @@ -211,7 +217,7 @@ i = self._k.index(key) k = self._k.pop(i) v = self._v.pop(i) - #exec('del self.%s')%(key) + #exec('del self.%s')%(key) return (k, v) def keys(self): @@ -226,8 +232,10 @@ def isempty(x): ''' - returns true if object is + -infinity, NaN, None, '', has length 0, or is an - array/matrix composed only of such components (includes mixtures of "empty" types)''' + returns true if object is +/-infinity, NaN, None, '', has length 0, or is + an array/matrix composed only of such components (includes mixtures of + "empty" types) + ''' if type(x) in [list, np.ndarray, tuple]: if np.size(x) == 0: @@ -261,8 +269,11 @@ def fieldnames(x, ignore_internals=True): - '''returns a list of fields of x - ignore_internals ignores all fieldnames starting with '_' and is True by default''' + ''' + returns a list of fields of x + ignore_internals ignores all fieldnames starting with '_' and is True by + default + ''' result = list(vars(x).keys()) if ignore_internals: @@ -272,8 +283,11 @@ def isfield(x, y, ignore_internals=True): - '''is y is a field of x? - ignore_internals ignores all fieldnames starting with '_' and is True by default''' + ''' + is y is a field of x? + ignore_internals ignores all fieldnames starting with '_' and is True by + default + ''' return str(y) in fieldnames(x, ignore_internals) @@ -280,7 +294,8 @@ def fileparts(x): ''' given: "path/path/.../file_name.ext" - returns: [path, file_name, ext] (list of strings)''' + returns: [path, file_name, ext] (list of strings) + ''' try: a = x[:x.rindex('/')] #path b = x[x.rindex('/') + 1:] #full filename @@ -296,17 +311,17 @@ def fullfile(*args): ''' - use: + usage: + fullfile(path, path, ... , file_name + ext) - fullfile(path, path, ... , file_name + ext) returns: "path/path/.../file_name.ext" with all arguments as strings with no "/"s regarding extensions and the '.': - as final arguments ('file.doc') or ('file' + '.doc') will work - ('final', '.doc'), and the like, will not (you'd get 'final/.doc') -''' + as final arguments ('file.doc') or ('file' + '.doc') will work + ('final', '.doc'), and the like, will not (you'd get 'final/.doc') + ''' result = str(args[0]) for i in range(len(args[1:])): # if last argument wasn't empty, add a '/' between it and the next argument @@ -320,8 +335,11 @@ def findline(fidi, s): ''' returns full first line containing s (as a string), or None + Note: will include any newlines or tabs that occur in that line, - use str(findline(f, s)).strip() to remove these, str() in case result is None''' + use str(findline(f, s)).strip() to remove these, str() in case result is + None + ''' for line in fidi: if s in line: return line @@ -330,20 +348,21 @@ def empty_nd_list(shape, filler=0., as_numpy_ndarray=False): ''' -returns a python list of the size/shape given (shape must be int or tuple) + returns a python list of the size/shape given (shape must be int or tuple) the list will be filled with the optional second argument filler is 0.0 by default - as_numpy_ndarray will return the result as a numpy.ndarray and is False by default + as_numpy_ndarray will return the result as a numpy.ndarray and is False by + default - Note: the filler must be either None/np.nan/float('NaN'), float/double, or int - other numpy and float values such as +/- np.inf will also work + Note: the filler must be either None/np.nan/float('NaN'), float/double, or + int. other numpy and float values such as +/- np.inf will also work -use: - empty_nd_list((5, 5), 0.0) # returns a 5x5 matrix of 0.0's - empty_nd_list(5, None) # returns a 5 long array of NaN -''' + Usage: + empty_nd_list((5, 5), 0.0) # returns a 5x5 matrix of 0.0's + empty_nd_list(5, None) # returns a 5 long array of NaN + ''' result = np.empty(shape) result.fill(filler) if not as_numpy_ndarray: Index: ../trunk-jpl/src/m/mech/analyticaldamage.py =================================================================== --- ../trunk-jpl/src/m/mech/analyticaldamage.py (revision 25064) +++ ../trunk-jpl/src/m/mech/analyticaldamage.py (revision 25065) @@ -1,6 +1,6 @@ import numpy as np + from averaging import averaging -#from plotmodel import plotmodel from thomasparams import thomasparams Index: ../trunk-jpl/src/m/plot/applyoptions.m =================================================================== --- ../trunk-jpl/src/m/plot/applyoptions.m (revision 25064) +++ ../trunk-jpl/src/m/plot/applyoptions.m (revision 25065) @@ -307,7 +307,6 @@ end end - %text (default value is empty, not NaN...) if exist(options,'text'); textstring=getfieldvalue(options,'text'); @@ -351,7 +350,7 @@ scaleruler(options); end -%streamliness +%streamlines if exist(options,'streamlines'), plot_streamlines(md,options); end Index: ../trunk-jpl/src/m/plot/plot_manager.py =================================================================== --- ../trunk-jpl/src/m/plot/plot_manager.py (revision 25064) +++ ../trunk-jpl/src/m/plot/plot_manager.py (revision 25065) @@ -1,14 +1,3 @@ - -from checkplotoptions import checkplotoptions -from plot_mesh import plot_mesh -from plot_BC import plot_BC -from plot_elementnumbering import plot_elementnumbering -from plot_vertexnumbering import plot_vertexnumbering -from processmesh import processmesh -from processdata import processdata -from plot_unit import plot_unit -from applyoptions import applyoptions - try: from osgeo import gdal overlaysupport = True @@ -16,8 +5,17 @@ print('osgeo/gdal for python not installed, overlay plots are not enabled') overlaysupport = False +from applyoptions import applyoptions +from checkplotoptions import checkplotoptions +from plot_BC import plot_BC +from plot_mesh import plot_mesh +from plot_elementnumbering import plot_elementnumbering if overlaysupport: from plot_overlay import plot_overlay +from plot_unit import plot_unit +from plot_vertexnumbering import plot_vertexnumbering +from processdata import processdata +from processmesh import processmesh def plot_manager(md, options, fig, axgrid, gridindex): @@ -26,9 +24,12 @@ 'fig' is a handle to the figure instance created by plotmodel. - 'ax' is a handle to the axes instance created by plotmodel. This is + 'axgrid' is a handle to the axes instance created by plotmodel. This is currently generated using matplotlib's AxesGrid toolkit. + 'gridindex' is passed through to specialized plot* functions and + applyoptions. + Usage: plot_manager(md, options, fig, ax) @@ -36,9 +37,14 @@ ''' #parse options and get a structure of options options = checkplotoptions(md, options) + #get data to be plotted data = options.getfieldvalue('data') + #add ticklabel has a default option + # + # TODO: Check why we are doing this and if it is absolutely necessary (see + # also src/m/plot/plot_mesh.py, src/m/plot/applyoptions.py) options.addfielddefault('ticklabels', 'on') ax = axgrid[gridindex] @@ -57,8 +63,7 @@ if isinstance(data, str): if data == 'mesh': plot_mesh(md, options, fig, axgrid, gridindex) - - #fig.delaxes(fig.axes[1]) # hack to remove colorbar after the fact + #fig.delaxes(fig.axes[1]) # hack to remove colorbar after the fact return elif data == 'BC': plot_BC(md, options, fig, axgrid, gridindex) Index: ../trunk-jpl/src/m/plot/plot_overlay.py =================================================================== --- ../trunk-jpl/src/m/plot/plot_overlay.py (revision 25064) +++ ../trunk-jpl/src/m/plot/plot_overlay.py (revision 25065) @@ -1,24 +1,31 @@ -import numpy as np -from processmesh import processmesh -from processdata import processdata -from xy2ll import xy2ll +import os + import matplotlib.pyplot as plt import matplotlib as mpl -import os try: from mpl_toolkits.basemap import Basemap except ImportError: print('Basemap toolkit not installed') +import numpy as np try: from osgeo import gdal except ImportError: print('osgeo/gdal for python not installed, plot_overlay is disabled') +from processdata import processdata +from processmesh import processmesh +from xy2ll import xy2ll + def plot_overlay(md, data, options, ax): ''' - Function for plotting a georeferenced image. This function is called - from within the plotmodel code. + Function for plotting a georeferenced image. Called from plot_manager by + call to plotmodel. + + Usage: + plot_overlay(md, data, options, ax) + + See also: PLOTMODEL ''' x, y, z, elements, is2d, isplanet = processmesh(md, [], options) @@ -80,7 +87,7 @@ plt.hist(arr.flatten(), bins=256, range=(0., 1.)) plt.title('histogram of overlay image, use for setting overlaylims') plt.show() - plt.sca(ax) # return to original axes / figure + plt.sca(ax) # return to original axes/figure # get parameters from cropped geotiff trans = gtif.GetGeoTransform() Index: ../trunk-jpl/src/m/plot/plotmodel.m =================================================================== --- ../trunk-jpl/src/m/plot/plotmodel.m (revision 25064) +++ ../trunk-jpl/src/m/plot/plotmodel.m (revision 25065) @@ -9,7 +9,7 @@ numberofplots=options.numberofplots; %get number of subplots -subplotwidth=ceil(sqrt(options.numberofplots)); +subplotwidth=ceil(sqrt(numberofplots)); %if nlines and ncols specified, then bypass. if ~exist(options.list{1},'nlines') & ~exist(options.list{1},'ncols') @@ -103,7 +103,7 @@ end end % }}} - %Use zbuffer renderer (snoother colors) and white background + %Use zbuffer renderer (smoother colors) and white background set(f,'Renderer','zbuffer','color',getfieldvalue(options.list{1},'figurebackgroundcolor','w')); %Go through all data plottable and close window if an error occurs @@ -124,5 +124,5 @@ rethrow(me); end else - error('plotmodel error message: no output data found. '); + error('plotmodel error message: no output data found.'); end Index: ../trunk-jpl/src/m/plot/plot_manager.m =================================================================== --- ../trunk-jpl/src/m/plot/plot_manager.m (revision 25064) +++ ../trunk-jpl/src/m/plot/plot_manager.m (revision 25065) @@ -14,9 +14,7 @@ %figure out if this is a special plot if ischar(data), - switch data, - case 'boundaries', plot_boundaries(md,options,subplotwidth,i); return; @@ -32,23 +30,18 @@ case 'highlightelements', plot_highlightelements(md,options,subplotwidth,i); return; - case 'qmumean', plot_qmumean(md,options,nlines,ncols,i); return; - case 'qmustddev', plot_qmustddev(md,options,nlines,ncols,i); return; - case 'qmuhistnorm', plot_qmuhistnorm(md,options,nlines,ncols,i); return; - case 'qmu_mass_flux_segments', plot_qmu_mass_flux_segments(md,options,nlines,ncols,i); return; - case 'part_hist', plot_parthist(md,options,nlines,ncols,i); return; @@ -120,10 +113,8 @@ case 'segments' plot_segments(md,options,subplotwidth,i,data) return - case 'quiver' data=[md.initialization.vx md.initialization.vy]; %Go ahead and try plot_unit - case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',... 'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',... 'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'}, @@ -137,13 +128,10 @@ return; case 'transient_results', plot_transient_results(md,options,subplotwidth,i); - case 'transient_field', plot_transient_field(md,options,subplotwidth,i); return; - otherwise, - if ismember(data,properties('model')), data=eval(['md.' data ';']); else @@ -158,19 +146,19 @@ return; end -%Figure out if this is a semi-transparent plot. +%Figure out if this is a Google Maps plot. if exist(options,'googlemaps'), plot_googlemaps(md,data,options,nlines,ncols,i); return; end -%Figure out if this is a semi-transparent plot. +%Figure out if this is a landsat plot. if getfieldvalue(options,'landsat',0), plot_landsat(md,data,options,nlines,ncols,i); return; end -%Figure out if this is a semi-transparent plot. +%Figure out if this is a gridded plot. if exist(options,'gridded'), plot_gridded(md,data,options,nlines,ncols,i); return; Index: ../trunk-jpl/src/m/shp/shpread.py =================================================================== --- ../trunk-jpl/src/m/shp/shpread.py (nonexistent) +++ ../trunk-jpl/src/m/shp/shpread.py (revision 25065) @@ -0,0 +1,136 @@ +from os import path +try: + import shapefile +except ImportError: + print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled") + +from pairoptions import pairoptions + + +def shpread(filename, *args): #{{{ + ''' + SHPREAD - read a shape file and build a dict + + This routine reads a shape file .shp and builds a dict array containing the + fields x and y corresponding to the coordinates, one for the filename + of the shp file, for the density, for the nodes, and a field closed to + indicate if the domain is closed. If this initial shapefile is point only, + the fields closed and points are ommited. + The first argument is the .shp file to be read and the second one + (optional) indicates if the last point shall be read (1 to read it, 0 not + to). + + The current implementation of shpread depends on PyShp. + + Usage: + dict = shpread(filename) + + Example: + From underling PyShp implementation, "The shapefile format is actually + a collection of three files. You specify the base filename of the + shapefile or the complete filename of any of the shapefile component + files."" + + dict = shpread('domainoutline.shp') + OR + dict = shpread('domainoutline.dbf') + OR + dict = shpread('domainoutline') + + "OR any of the other 5+ formats which are potentially part of a + shapefile. The library does not care about file extensions". We do, + however, check that a file with the base filename or base filename with + .shp extension exists. + + See also EXPDOC, EXPWRITEASVERTICES + + Sources: + - https://github.com/GeospatialPython/pyshp + ''' + + #recover options + options = pairoptions(*args) + + #some checks + if not (path.isfile(filename) or path.isfile(filename + '.shp')): + raise RuntimeError('shpread error message: file {} or {}.shp not found!'.format(filename, filename)) + + #read shapefile + sf = shapefile.Reader(filename) + + Dicts = [] + + #retrieve ID (if it exists) + fields = sf.fields + for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag" + if fields[i][0] == 'id': + name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag" + break + + shapes = sf.shapes() + for i in range(len(shapes)): + Dict = {} + shape = shapes[i] + if shape.shapeType == shapefile.POINT: + Dict['x'] = shape.points[0][0] + Dict['y'] = shape.points[0][1] + Dict['density'] = 1 + Dict['Geometry'] = 'Point' + elif shape.shapeType == shapefile.POLYLINE: + num_points = len(shape.points) + x = [] + y = [] + for j in range(num_points): + point = shape.points[j] + x.append(point[0]) + y.append(point[1]) + Dict['x'] = x + Dict['y'] = y + Dict['nods'] = num_points + Dict['density'] = 1 + Dict['closed'] = 1 + Dict['BoundingBox'] = shape.bbox + Dict['Geometry'] = 'Line' + elif shape.shapeType == shapefile.POLYGON: + num_points = len(shape.points) + x = [] + y = [] + for j in range(num_points): + point = shape.points[j] + x.append(point[0]) + y.append(point[1]) + Dict['x'] = x + Dict['y'] = y + Dict['nods'] = num_points + Dict['density'] = 1 + Dict['closed'] = 1 + Dict['BoundingBox'] = shape.bbox + Dict['Geometry'] = 'Polygon' + else: + # NOTE: We could do this once before looping over shapes as all + # shapes in the file must be of the same type, but we would + # need to have a second check anyway in order to know how to + # parse the points. So, let's just assume the file is not + # malformed. + # + raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName)) + name = '' + fields = sf.fields + for j in range(1, len(fields)): # skip over first field, which is "DeletionFlag" + fieldname = fields[j][0] + # 'id' field gets special treatment + if fieldname == 'id': + name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag" + else: + Dict[fieldname] = sf.record(j - 1)[0] + Dict['name'] = name + Dicts.append(Dict) + + invert = options.getfieldvalue('invert', 0) + if invert: + for i in range(len(Dicts)): + Dicts[i].x = np.flipud(Dicts[i].x) + Dicts[i].y = np.flipud(Dicts[i].y) + + return Dicts +#}}} Index: ../trunk-jpl/src/m/shp/shpread.m =================================================================== --- ../trunk-jpl/src/m/shp/shpread.m (revision 25064) +++ ../trunk-jpl/src/m/shp/shpread.m (revision 25065) @@ -1,22 +1,22 @@ function Struct=shpread(filename,varargin) -%SHPREAD - read a shape file and build a Structure +%SHPREAD - read a shape file and build a struct % -% This routine reads a shape file .shp and builds a Structure containing the -% fields x and y corresponding to the coordinates, one for the filename of -% the shp file, for the density, for the nodes, and a field closed to -% indicate if the domain is closed. -% If this initial shapefile is point only, the fields closed and -% points are ommited -% The first argument is the .shp file to be read and the second one (optional) -% indicates if the last point shall be read (1 to read it, 0 not to). +% This routine reads a shape file .shp and builds a structure array +% containing the fields x and y corresponding to the coordinates, one for the +% filename of the shp file, for the density, for the nodes, and a field +% closed to indicate if the domain is closed. If this initial shapefile is +% point only, the fields closed and points are omitted. +% The first argument is the .shp file to be read and the second one +% (optional) indicates if the last point shall be read (1 to read it, 0 not +% to). % -% Usage: -% Struct=shpread(filename) +% Usage: +% Struct=shpread(filename) % -% Example: -% Struct=shpread('domainoutline.shp') +% Example: +% Struct=shpread('domainoutline.shp') % -% See also EXPDOC, EXPWRITEASVERTICES +% See also EXPDOC, EXPWRITEASVERTICES %recover options options=pairoptions(varargin{:}); @@ -26,9 +26,6 @@ error(['shpread error message: file ' filename ' not found!']); end -%initialize number of profile -count=0; - %read shapefile shp=shaperead(filename); @@ -35,7 +32,7 @@ Struct=struct([]); fields=fieldnames(shp); for i=1:length(shp), - if strcmpi(shp(i).Geometry,'Polygon'), + if strcmpi(shp(i).Geometry,'Point'), x=shp(i).X'; y=shp(i).Y'; ids=find(isnan(x)); x(ids)=[]; y(ids)=[]; @@ -42,9 +39,7 @@ Struct(end+1).x=x; Struct(end).y=y; - Struct(end).nods=length(x); Struct(end).density=1; - Struct(end).closed=1; if isfield(shp,'id'), Struct(end).name=num2str(shp(i).id); else @@ -56,9 +51,7 @@ Struct(end).(field)=shp(i).(field); end end - end - - if strcmpi(shp(i).Geometry,'Line'), + elseif strcmpi(shp(i).Geometry,'Line'), x=shp(i).X'; y=shp(i).Y'; ids=find(isnan(x)); x(ids)=[]; y(ids)=[]; @@ -79,10 +72,7 @@ Struct(end).(field)=shp(i).(field); end end - end - - - if strcmpi(shp(i).Geometry,'Point'), + elseif strcmpi(shp(i).Geometry,'Polygon'), x=shp(i).X'; y=shp(i).Y'; ids=find(isnan(x)); x(ids)=[]; y(ids)=[]; @@ -89,7 +79,9 @@ Struct(end+1).x=x; Struct(end).y=y; + Struct(end).nods=length(x); Struct(end).density=1; + Struct(end).closed=1; if isfield(shp,'id'), Struct(end).name=num2str(shp(i).id); else Index: ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py =================================================================== --- ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 25064) +++ ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 25065) @@ -1,4 +1,3 @@ -import os import subprocess import numpy as np @@ -9,7 +8,7 @@ from pairoptions import * -def gmshplanet(*varargin): +def gmshplanet(*args): ''' GMSHPLANET - mesh generation for a sphere. Very specific code for gmsh from $ISSM_DIR/src/demos/simple_geo/sphere.geo @@ -33,7 +32,7 @@ raise RuntimeError('gmshplanet: Gmsh major version %s not supported!' % gmshmajorversion) #process options - options = pairoptions(*varargin) + options = pairoptions(*args) #options = deleteduplicates(options, 1) #recover parameters: Index: ../trunk-jpl/src/m/coordsystems/gdaltransform.py =================================================================== --- ../trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25064) +++ ../trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25065) @@ -5,31 +5,35 @@ from loadvars import * + def gdaltransform(x, y, proj_in, proj_out): #{{{ ''' GDALTRANSFORM - switch from one projection system to another - Usage: - [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out) + Usage: + [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out) - Example: - [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031') + Example: + [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031') - For reference: - EPSG: 4326 (lat, long) - EPSG: 3341 (Greenland, UPS 45W, 70N) - EPSG: 3031 (Antarctica, UPS 0E, 71S) + For reference: + EPSG: 4326 (lat, long) + EPSG: 3341 (Greenland, UPS 45W, 70N) + EPSG: 3031 (Antarctica, UPS 0E, 71S) - ll2xy default projection Antarctica: - +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 - ll2xy default projection Greenland: - +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 - Bamber's Greenland projection - +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 + ll2xy default projection Antarctica: + +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 + ll2xy default projection Greenland: + +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 + Bamber's Greenland projection + +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 - To get proj.4 string form EPSG, use gdalsrsinfo. Example: + To get proj.4 string form EPSG, use gdalsrsinfo. Example: - gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" + gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" + + TODO: + - Remove shlex and follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py ''' # Give ourselves unique file names @@ -41,7 +45,7 @@ file_in.write(b'%8g %8g\n' % (x.flatten(1) y.flatten(1))) file_in.close() - args = shlex.split('gdaltransform -s_srs %s -t_srs %s < %s > %s' % (proj_in, proj_out, filename_in, filename_out)) + args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out)) subprocess.check_call(args) # Will raise CalledProcessError if return code is not 0 A = loadvars(filename_out) @@ -52,4 +56,6 @@ os.remove(filename_in) os.remove(filename_out) + + return [xout, yout] #}}} Index: ../trunk-jpl/src/m/geometry/polyarea.py =================================================================== --- ../trunk-jpl/src/m/geometry/polyarea.py (nonexistent) +++ ../trunk-jpl/src/m/geometry/polyarea.py (revision 25065) @@ -0,0 +1,28 @@ +import math + +import numpy as np + + +def polyarea: #{{{ + ''' + POLYAREA - returns the area of the 2-D polygon defined by the vertices in + lists x and y + + Partial implementation of MATLAB's polyarea function. If x and y are + lists of the same length, then polyarea returns the scalar area of the + polygon defined by x and y. + + Usage: + a = polyarea(x, y) + + Sources: + - https://www.mathworks.com/help/matlab/ref/polyarea.html + - https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates + - https://en.wikipedia.org/wiki/Shoelace_formula + + TODO: + - Test that output falls within some tolerance of MATLAB's polyarea + function. + ''' + return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) +#}}} Index: ../trunk-jpl/src/m/classes/boundary.py =================================================================== --- ../trunk-jpl/src/m/classes/boundary.py (nonexistent) +++ ../trunk-jpl/src/m/classes/boundary.py (revision 25065) @@ -0,0 +1,176 @@ +import math + +import numpy as np + +from epsg2proj import epsg2proj +from pairoptions import pairoptions +from shpread import shpread + + +class boundary(object): #{{{ + ''' + BOUNDARY class definition + + Usage: + boundary = boundary() + ''' + + def __init__(self, *args): #{{{ + self.boundaries = [] + self.name = '' + self.continent = '' + self.proj = epsg2proj(4326) + + self.setdefaultparameters() + + if len(args): + options = pairoptions(*args) + + #recover field values + self.shppath = options.getfieldvalue('shppath', '') + self.shpfilename = options.getfieldvalue('shpfilename', '') + self.orientation = options.getfieldvalue('orientation', 'normal') + + #figure out projection string: + if options.exist('epsg'): + if options.exist('proj'): + raise RuntimeError('Error in constructor for boundary {}. Cannot supply epsg and proj at the same time!'.format(self.shppath)) + epsg = options.getfieldvalue('epsg') + proj = epsg2proj(epsg) # convert to PROJ.4 format + elif options.exist('proj'): + if options.exist('epsg'): + raise RuntimeError('Error in constructor for boundary {}. Cannot supply proj and epsg at the same time!'.format(self.shppath)) + proj = options.getfieldvalue('proj') + else: + proj = '+proj=longlat +datum=WGS84 +no_defs' + + self.proj = proj + #}}} + + def __repr__(self): #{{{ + s = ' boundary parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'shppath', 'path to filename for this boundary')) + s += '{}\n'.format(fielddisplay(self, 'shpfilename', 'shape file name')) + s += '{}\n'.format(fielddisplay(self, 'orientation', 'orientation (default is \'normal\', can be \'reverse\')')) + s += '{}\n'.format(fielddisplay(self, 'proj', 'shape file projection string (follows PROJ.4 standard)')) + + return s + #}}} + + def setdefaultparameters(self): #{{{ + self.shppath = '' + self.shpfilename = '' + self.orientation = 'normal' + self.proj = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326 + + return self + #}}} + + def name(self): #{{{ + output = self.shpfilename + + return output + #}}} + + def edges(self): #{{{ + #read domain + path, name, ext = fileparts(shp.shpfilename) + if ext != 'shp': + ext = 'shp' + output = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + + #do we reverse? + if self.orientation == 'reverse': + output.x = np.flipud(output.x) + output.y = np.flipud(output.y) + #}}} + + def plot(self, *args): #{{{ + #recover options + options = pairoptions(*args) + + #parse input + # TODO: Sort out what options we should set (see + # src/m/classes/boundary.m) + + if options.exist('epsg'): + if options.exist('proj'): + raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection') + proj = epsg2proj(options.getfieldvalue('epsg')) + elif options.exist('proj'): + proj = options.getfieldvalue('proj') + else: + proj = epsg2proj(4326) + + #read domain + path, name, ext = fileparts(shp.shpfilename) + if ext != 'shp': + ext = 'shp' + domain = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + + #convert boundary to another reference frame #{{{ + for i in range(len(domain)): + try: + x, y = gdaltransform(domain[i].x, domain[i].y, self.proj, proj) + except error as e: + print(e) + print(self) + + # TODO: Figure out how to recover figure here: do we pass 'fig' and + # 'ax' in args? + #for i in range(len(domain)): + # x = domain[i].x * unitmultiplier + # y = domain[i].y * unitmultiplier + # if len(x) == 1: + #}}} + + #TODO: Finish translating from MATLAB after test2004.py runs without plot + #}}} + + def checkconsistency(self, *args): #{{{ + #recover options + options = pairoptions(*args) + tolerance = getfieldvalue(options, 'tolerance', 1e-5) + + #recover edges + edges = self.edges() + + if edges.Geometry == 'Point': + return + else: + #check that we don't have two identical vertices + x = edges.x + y = edges.y + distance = math.sqrt((x[:-2] - x[1:-1]) ** 2 + (y[:-2] - y[1:-1]) ** 2) + dmax = distance.max() + isidentical = np.where(np.asarray(distance) < dmax * tolerance) + for elem in isidentical: # distance should be a 1D array, but return from np.where is a tuple of arrays + if len(elem) != 0: + raise Exception('boundary {} has two vertices extermely close to one another'.format(shpfilename)) + + def plot3d(self, *args): #{{{ + #recover options + options = pairoptions(*args) + + #parse input + # TODO: Sort out what options we should set (see + # src/m/classes/boundary.m) + + if options.exist('epsg'): + if options.exist('proj'): + raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection') + proj = epsg2proj(options.getfieldvalue('epsg')) + elif options.exist('proj'): + proj = options.getfieldvalue('proj') + else: + proj = epsg2proj(4326) + + #read domain + path, name, ext = fileparts(shp.shpfilename) + if ext != 'shp': + ext = 'shp' + domain = shpread('{}/{}.{}'.format(self.shppath, name, ext)) + + #TODO: Finish translating from MATLAB after test2004.py runs without plot + #}}} +#}}} Index: ../trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 25064) +++ ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 25065) @@ -162,10 +162,10 @@ list=unique(list); end % }}} function addbasin(self,bas) % {{{ - if ~strcmpi(class(bas),'basin') - error('addbasin method only takes a ''basin'' class object as input'); - end; - self.basins{end+1}=bas; + if ~strcmpi(class(bas),'basin') + error('addbasin method only takes a ''basin'' class object as input'); + end; + self.basins{end+1}=bas; end % }}} function intersections(self,varargin) % {{{ Index: ../trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 25064) +++ ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 25065) @@ -10,21 +10,22 @@ from private import private from TwoDToThreeD import TwoDToThreeD + class sealevelmodel(object): ''' SEALEVELMODEL class definition - Usage: - slm = sealevelmodel(*args) + Usage: + slm = sealevelmodel(*args) - where args is a variable list of options + where args is a variable list of options - Example: - slm = sealevelmodel( - 'icecap', md_greenland, - 'icecap', md_antarctica, - 'earth', md_earth - ) + Example: + slm = sealevelmodel( + 'icecap', md_greenland, + 'icecap', md_antarctica, + 'earth', md_earth + ) ''' def __init__(self, *args): #{{{ @@ -58,11 +59,11 @@ #}}} def __repr__(self): # {{{ - s = '%s\n' % fielddisplay(self, 'icecaps', 'ice caps') - s += '%s\n' % fielddisplay(self, 'earth', 'earth') - s += '%s\n' % fielddisplay(self, 'settings', 'settings properties') - s += '%s\n' % fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...') - s += '%s\n' % fielddisplay(self, 'miscellaneous', 'miscellaneous fields') + s = '{}\n'.format(fielddisplay(self, 'icecaps', 'ice caps')) + s += '{}\n'.format(fielddisplay(self, 'earth', 'earth')) + s += '{}\n'.format(fielddisplay(self, 'settings', 'settings properties')) + s += '{}\n'.format(fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...')) + s += '{}\n'.format(fielddisplay(self, 'miscellaneous', 'miscellaneous fields')) #}}} def setdefaultparameters(self): # {{{ @@ -84,7 +85,7 @@ # Is the coupler turned on? for i in range(len(slm.icecaps)): if slm.icecaps[i].transient.iscoupler == 0: - warnings.warn('sealevelmodel checkconsistency error: icecap model %s should have the transient coupler option turned on!' % slm.icecaps[i].miscellaneous.name) + warnings.warn('sealevelmodel checkconsistency error: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name)) if slm.earth.transient.iscoupler == 0: warnings.warn('sealevelmodel checkconsistency error: earth model should have the transient coupler option turned on!') @@ -92,18 +93,18 @@ # Check that the transition vectors have the right size for i in range(len(slm.icecaps)): if slm.icecaps[i].mesh.numberofvertices != len(slm.earth.slr.transitions[i]): - raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: %i name: %s' % (i, slm.icecaps[i].miscellaneous.name)) + raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name)) # Check that run frequency is the same everywhere for i in range(len(slm.icecaps)): if slm.icecaps[i].slr.geodetic_run_frequency != slm.earth.geodetic_run_frequency: - raise RuntimeError('sealevelmodel checkconsistency error: icecap model %s should have the same run frequency as earth!' % slm.icecaps[i].miscellaneous.name) + raise RuntimeError('sealevelmodel checkconsistency error: icecap model {} should have the same run frequency as earth!'.format(slm.icecaps[i].miscellaneous.name)) # Make sure steric_rate is the same everywhere for i in range(len(slm.icecaps)): md = slm.icecaps[i] if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []: - raise RuntimeError('steric rate on ice cap %s is not the same as for the earth' % md.miscellaneous.name) + raise RuntimeError('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name)) #}}} def mergeresults(self): # {{{ @@ -137,7 +138,7 @@ def listcaps(self): # {{{ for i in range(len(self.icecaps)): - print('%i: %s' % (i, self.icecaps[i].miscellaneous.name)) + print('{}: {}'.format(i, self.icecaps[i].miscellaneous.name)) #}}} def continents(self): # {{{ @@ -184,7 +185,7 @@ yei = mdi.mesh.y[mdi.mesh.elements] * onesmatrix / 3 zei = mdi.mesh.z[mdi.mesh.elements] * onesmatrix / 3 - print('Computing vertex intersections for basin %s' % self.basins[i].name) + print('Computing vertex intersections for basin {}'.format(self.basins[i].name)) 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)) self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force)) Index: ../trunk-jpl/src/m/classes/basin.py =================================================================== --- ../trunk-jpl/src/m/classes/basin.py (nonexistent) +++ ../trunk-jpl/src/m/classes/basin.py (revision 25065) @@ -0,0 +1,214 @@ +import math + +from boundary import boundary +from epsg2proj import epsg2proj +from helpers import fileparts +from laea import laea +from pairoptions import pairoptions +from polyarea import polyarea +from shpread import shpread + + +class basin(object): #{{{ + ''' + BASIN class definition + + Usage: + basin = basin() + ''' + + def __init__(self, *args): #{{{ + self.boundaries = [] + self.name = '' + self.continent = '' + self.proj = epsg2proj(4326) + + self.setdefaultparameters() + + if len(args): + options = pairoptions(*args) + + #recover field values + self.boundaries = options.getfieldvalue('boundaries', []) + self.name = options.getfieldvalue('name', '') + self.continent = options.getfieldvalue('continent', '') + + # figure out projection string + if options.exist('epsg'): + if options.exist('proj'): + raise RuntimeError('Error in constructor for basin %s. Cannot supply epsg and proj at the same time!'.format(self.name)) + epsg = options.getfieldvalue('epsg', '') + proj = epsg2proj(epsg)#convert to PROJ.4 format + elif options.exist('proj'): + if options.exist('epsg'): + raise RuntimeError('Error in constructor for basin {}. Cannot supply proj and epsg at the same time!'.format(self.name)) + proj = options.getfieldvalue('proj', '') + else: + proj = '+proj=longlat +datum=WGS84 +no_defs' + + self.proj = proj + #}}} + + def __repr__(self): # {{{ + s = ' basin parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'continent', 'continent name')) + s += '{}\n'.format(fielddisplay(self, 'name', 'basin name')) + s += '{}\n'.format(fielddisplay(self, 'proj', 'basin projection string (follows PROJ.4 standard')) + s += '{}\n'.format(fielddisplay(self, 'boundaries','list of boundary objects')) + for i in range(len(self.boundaries)): + s += ' boundary #{}: {}\n'.format((i + 1), self.boundaries[i]) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + self.name = '' + self.continent = '' + self.proj = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326 + self.boundaries = [] + + return self + #}}} + + def isnameany(self, *args): #{{{ + for arg in args: + if arg == self.name: + return 1 + return 0 + #}}} + + def iscontinentany(self, *args): #{{{ + for arg in args: + if arg == self.continent: + return 1 + return 0 + #}}} + + def outputname(self, *args): #{{{ + #recover options + options = pairoptions(*args) + extension = options.getfieldvalue('extension', 1) + + path, name, ext = fileparts(*args) + if extension: + output = '{}{}'.format(name, ext) + else: + output = name + + return output + #}}} + + def plot(self, *args): #{{{ + #add option + for i in range(len(self.boundaries)): + self.boundaries[i].plot('proj', self.proj, *args) + #}}} + + def plot3d(self, *args): #{{{ + #add option + for i in range(len(self.boundaries)): + self.boundaries[i].plot3d(*args) + #}}} + + def contour(self, *args): #{{{ + #recover options + options = pairoptions(*args) + x = [] + y = [] + + #go through boundaries, recover edges and project them in the basin epsg reference frame + for i in range(len(self.boundaries)): + boundary = self.boundaries[i] + contour = boundary.edges() + contour.x, contour.y = gdaltransform(contour.x, contour.y, boundary.proj, self.proj) + x = [[x], [contour.x]] + y = [[y], [contour.y]] + + #close the contour + if x[-1] != x[0] or y[-1] != y[0]: + x.append(x[0]) + y.append(y[0]) + + setattr(out, 'x', x) + setattr(out, 'y', y) + setattr(out, 'nods', len(x)) + + return out + #}}} + + def checkconsistency(self, *args): #{{{ + #recover options + options = pairoptions(*args) + + #figure out if two boundaries are identical + for i in range(len(self.boundaries)): + namei = self.boundaries[i].shpfilename + for j in range(1, len(self.boundaries)): + namej = self.boundaries[j].shpfilename + if namei == namej: + raise RuntimeError('Basin {} has two identical boundaries named {}'.format(self.name, namei)) + + #go through boundaries and check their consistency + for i in range(len(self.boundaries)): + boundary == self.boundaries[i] + boundary.checkconsistency() + #}}} + + def contourplot(self, *args): #{{{ + contour = self.contour() + plot(contour.x, contour.y, 'r*-') + #}}} + + def shapefilecrop(self, *args): #{{{ + #recover options + options = pairoptions(*args) + threshold = options.getfieldvalue('threshold', .65) #.65 degrees lat, long + inshapefile = options.getfieldvalue('shapefile') + outshapefile = options.getfieldvalue('outshapefile', '') + + if options.exist('epsgshapefile') and options.exist('projshapefile'): + raise RuntimeError('Basin shapefilecrop error messgae: cannot specify epsgshapefile and projshapefile at the same time') + + if options.exist('epsgshapefile'): + projshapefile = epsg2proj(options.getfieldvalue('epsgshapefile')) + elif options.exist('projshapefile'): + projshapefile = options.getfieldvalue('projshapefile') + else: + raise RuntimeError('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified') + + #create list of contours that have critical length > threshold (in lat, long) + contours = shpread(inshapefile) + llist = [] + for i in range(len(contours)): + contour = contours[i] + carea = polyarea(contour.x, contour.y) + clength = math.sqrt(carea) + if clength < threshold: + llist = np.concatenate(llist, i) + setattr(contours, llist, []) + + #project onto reference frame + for i in range(len(contours)): + h = contours[i] + h.x, h.y = gdaltransform(h.x, h.y, projshapefile, self.proj) + contours[i].x = h.x + contours[i].y = h.y + + #only keep the contours that are inside the basin (take centroids) + ba = self.contour() + flags = np.zeros(len(contours)) + for i in range(len(contours)) + h = contours[i] + isin = inpolygon(h.x, h.y, ba.x, ba.y) + if len(np.where(isin == 0, isin, isin)): + flags[i] = 1 + pos = flags.nonzero() + contours[pos] = [] + + #Two options + if outputshapefile == '': + output = contours + else: + shpwrite(contours, outputshapefile) + #}}} +#}}} Index: ../trunk-jpl/src/m/classes/boundary.m =================================================================== --- ../trunk-jpl/src/m/classes/boundary.m (revision 25064) +++ ../trunk-jpl/src/m/classes/boundary.m (revision 25065) @@ -19,35 +19,32 @@ end methods function self = boundary(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - otherwise - self=setdefaultparameters(self); - options=pairoptions(varargin{:}); - - %recover field values: - self.shppath=getfieldvalue(options,'shppath',''); - self.shpfilename=getfieldvalue(options,'shpfilename',''); - self.orientation=getfieldvalue(options,'orientation','normal'); + self=setdefaultparameters(self); - %figure out projection string: + if nargin > 0, + options=pairoptions(varargin{:}); + + %recover field values: + self.shppath=getfieldvalue(options,'shppath',''); + self.shpfilename=getfieldvalue(options,'shpfilename',''); + self.orientation=getfieldvalue(options,'orientation','normal'); + + %figure out projection string: + if exist(options,'epsg'), + if exist(options,'proj'), + error(['Error in constructor for boundary ' self.shppath '. Cannot supply epsg and proj at the same time!']); + end + epsg=getfieldvalue(options,'epsg'); + proj=epsg2proj(epsg); % convert to PROJ.4 format + elseif exist(options,'proj'), if exist(options,'epsg'), - if exist(options,'proj'), - error(['Error in constructor for boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']); - end - epsg=getfieldvalue(options,'epsg'); - %convert into proj4: - proj=epsg2proj(epsg); - elseif exist(options,'proj'), - if exist(options,'epsg'), - error(['Error in constructor for boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']); - end - proj=getfieldvalue(options,'proj'); - else - proj= '+proj=longlat +datum=WGS84 +no_defs'; + error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']); end - self.proj=proj; + proj=getfieldvalue(options,'proj'); + else + proj= '+proj=longlat +datum=WGS84 +no_defs'; + end + self.proj=proj; end end % }}} function self = setdefaultparameters(self) % {{{ @@ -54,7 +51,7 @@ self.shppath=''; self.shpfilename=''; self.orientation='normal'; - self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326 + self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326 end % }}} function disp(self) % {{{ disp(sprintf(' boundary parameters:')); @@ -69,19 +66,19 @@ output=self.shpfilename; end % }}} function output=edges(self) % {{{ - - %read domain: - [path,name,ext]=fileparts(self.shpfilename); - if ~strcmpi(ext,'shp'), - ext='shp'; - end - output=shpread([self.shppath '/' name '.' ext]); - %do we reverse? - if strcmpi(self.orientation,'reverse'), - output.x=flipud(output.x); - output.y=flipud(output.y); - end + %read domain: + [path,name,ext]=fileparts(self.shpfilename); + if ~strcmpi(ext,'shp'), + ext='shp'; + end + output=shpread([self.shppath '/' name '.' ext]); + + %do we reverse? + if strcmpi(self.orientation,'reverse'), + output.x=flipud(output.x); + output.y=flipud(output.y); + end end % }}} function plot(self,varargin) % {{{ %recover options @@ -100,10 +97,9 @@ fontsize=getfieldvalue(options,'fontsize',10); label=getfieldvalue(options,'label','on'); - if (exist(options,'epsg') & exist(options,'proj')), - error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection'); - end - if exist(options,'epsg'), + if (exist(options,'epsg'), + if exist(options,'proj')), + error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection'); proj=epsg2proj(getfieldvalue(options,'epsg')); elseif exist(options,'proj'), proj=getfieldvalue(options,'proj'); @@ -118,7 +114,7 @@ end domain=shpread([self.shppath '/' name '.' ext]); - %convert boundary to another reference frame: {{{ + %convert boundary to another reference frame: {{{ for i=1:length(domain), try, [x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj); @@ -168,7 +164,7 @@ dmax=max(distance); isidentical=find(distance