Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25035) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25035) @@ -0,0 +1,191 @@ +import warnings + +import numpy as np + +from generic import generic +from issmsettings import issmsettings +from meshintersect3d import meshintersect3d +from miscellaneous import miscellaneous +from pairoptions import pairoptions +from private import private +from TwoDToThreeD import TwoDToThreeD + +class sealevelmodel(object): + ''' + SEALEVELMODEL class definition + + Usage: + slm = sealevelmodel(*args) + + where args is a variable list of options + + Example: + slm = sealevelmodel( + 'icecap', md_greenland, + 'icecap', md_antarctica, + 'earth', md_earth + ) + ''' + + def __init__(self, *args): #{{{ + self.icecaps = [] # list of land/ice models; name should be changed later + self.earth = 0 # model for the whole earth + self.basins = [] # list of basins, matching icecaps, where shapefile info is held + self.cluster = 0 + self.miscellaneous = 0 + self.settings = 0 + self.private = 0 + self.mergedcaps = 0 + self.transitions = [] + self.eltransitions = [] + self.planet = '' + + # Create a default object + self.setdefaultparameters() + + if len(args): + # Use provided options to set fields + options = pairoptions(*args) + + # Recover all the icecap models + self.icecaps = options.getfieldvalue('ice_cap', []) + + # Recover the earth models + self.earth = options.getfieldvalue('earth', 0) + + # Set planet type + self.planet = options.getfieldvalue('planet', 'earth') + #}}} + + 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') + #}}} + + def setdefaultparameters(self): # {{{ + self.icecaps = [] + self.earth = [] + self.cluster = generic() + self.miscellaneous = miscellaneous() + self.settings = issmsettings() + self.private = private() + self.transitions = [] + self.eltransitions = [] + self.planet = 'earth' + + return self + #}}} + + @staticmethod + def checkconsistency(slm, solutiontype): # {{{ + # 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) + + if slm.earth.transient.iscoupler == 0: + warnings.warn('sealevelmodel checkconsistency error: earth model should have the transient coupler option turned on!') + + # 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)) + + # 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) + + # 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) + #}}} + + def mergeresults(self): # {{{ + champs = fieldnames(self.icecaps[i].results.TransientSolution) + for i in range(len(self.mergedcaps / 2)): + md = self.mergedcaps[2 * i] + trans = self.mergedcaps[2 * i + 1] + for j in range(len(self.icecaps[0].results.TransientSolution)): + for k in range(len(champs)): + if isinstance(getattr(icecaps[0].results.TransientSolution[j], champs[k]), float): + # Vertex or element? + if len(getattr(icecaps[0].results.TransientSolution[j], champs[k]) == icecaps[0].mesh.numberofvertices): + setattr(md.results.TransientSolution[j], champs[k], np.zeros(md.mesh.numberofvertices)) + for l in range(len(trans)): + resultcap = getattr(icecaps[l].results.TransientSolution[j], champs[k]) + setattr(getattr(md.results.TransientSolution[j], champs[k]), trans[l], resultcap) + else: + if champs[k] == 'IceVolume' or champs[k] == 'IceVolumeAboveFlotation': + setattr(md.results.TransientSolution, champs[k], 0) + for l in range(len(trans)): + resultcap = getattr(icecaps[l].results.TransientSolution[j], champs[k]) + setattr(md.results.TransientSolution[j], champs[k], getattr(md.results.TransientSolution[j], champs[k]) + resultcap) + elif champs[k] == 'time': + setattr(md.results.TransientSolution[j], champs[k], getattr(icecaps[0].results.TransientSolution[j], champs[k])) + else: + continue + else: + continue + self.mergedcaps[2 * i] = md + #}}} + + def listcaps(self): # {{{ + for i in range(len(self.icecaps)): + print('%i: %s' % (i, self.icecaps[i].miscellaneous.name)) + #}}} + + def continents(self): # {{{ + list = [] + for i in range(len(self.basins)): + list.append = self.basins[i].continent + return np.unique(list) + #}}} + + def basinsfromcontinent(self, continent): # {{{ + list = [] + for i in range(len(self.icecaps)): + if self.basins[i].continent == continent: + list.append = self.basins[i].name + return np.unique(list) + #}}} + + def addbasin(self, bas): # {{{ + if bas.__class__.__name__ != 'basin': + raise RuntimeError('addbasin method only takes a \'basin\' class object as input') + self.basins.append(bas) + #}}} + + def intersections(self, *args): #{{{ + options = pairoptions(*args) + force = options.getfieldvalue('force', 0) + + # Initialize, to avoid issues of having more transitions than meshes + self.transitions = [] + self.eltransitions = [] + + # For elements + onesmatrix = np.array([[1], [1], [1]]) + xe = self.earth.mesh.x[self.earth.mesh.elements] * onesmatrix / 3 + ye = self.earth.mesh.y[self.earth.mesh.elements] * onesmatrix / 3 + ze = self.earth.mesh.z[self.earth.mesh.elements] * onesmatrix / 3 + + for i in range(len(self.icecaps)): + mdi = self.icecaps[i] + mdi = TwoDToThreeD(mdi, self.planet) + + # For elements + zei = mdi.mesh.x[mdi.mesh.elements] * onesmatrix / 3 + 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) + + 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: /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m (revision 25034) +++ /issm/trunk-jpl/src/m/coordsystems/epsg2proj.m (revision 25035) @@ -1,14 +1,15 @@ function string = epsg2proj(epsg) -%use gdalsrsinfo to provide Proj.4 compatible string from EPSG code +%FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string from +%EPSG code % % Usage: % proj4string = epsg2proj(4326); % -% Example: -% proj4string = epsg2proj(4326); -% return proj4string='+proj=longlat +datum=wgs84 +no_defs ' +% Example: +% proj4string = epsg2proj(4326); +% return proj4string='+proj=longlat +datum=wgs84 +no_defs' % - [s,r]=system(['gdalsrsinfo epsg:' num2str(epsg) ' | grep PROJ.4 | sed "s/PROJ.4 : //"']); + [s,r]=system(['gdalsrsinfo epsg:' num2str(epsg) ' | grep PROJ.4 | sed "s/PROJ.4 : //"']); if s~=0, Index: /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py (revision 25035) +++ /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py (revision 25035) @@ -0,0 +1,29 @@ +import shlex +import subprocess + +def epsg2proj(epsg): #{{{ + ''' + FUNCTION EPSG2PROJ - uses gdalsrsinfo to provide PROJ.4 compatible string + from EPSG code + + Usage: + proj4string = epsg2proj(4326) + + Example: + proj4string = epsg2proj(4326) + return proj4string = '+proj=longlat +datum=wgs84 +no_defs' + + TODO: + - Implement try/catch for proc.communicate() + - In case of Python 2, except socket.timeout: https://docs.python.org/3/library/socket.html?highlight=socket%20timeout#socket.timeout + - In case of Python 3, except TimeoutExpired: https://docs.python.org/3/library/subprocess.html#subprocess.SubprocessError + ''' + args = 'gdalsrsinfo epsg:%s | grep PROJ.4 | sed "s/PROJ.4 : //"' % epsg + proc = subprocess.Popen(args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT) + outs, errs = proc.communicate() + + if errs != None: + raise RuntimeError('epsg2proj: call to gdalsrsinfo failed: %s' % errs) + + return outs +#}}} Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25034) +++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m (revision 25035) @@ -2,23 +2,25 @@ %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] = gdaltransform(md.mesh.long,md.mesh.lat,'EPSG:4326','EPSG:3031') +% Example: +% [x,y] = gdaltransform(md.mesh.long,md.mesh.lat,'EPSG:4326','EPSG:3031') % -% For reference: -% EPSG: 4326 (lat,long) -% EPSG: 3411 (Greenland, UPS 45W, 70N) -% EPSG: 3031 (Antarctica, UPS 0E, 71S) +% For reference: +% EPSG: 4326 (lat,long) +% EPSG: 3411 (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 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 +% 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 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: gdalsrsinfo "EPSG:4326" | grep "PROJ.4" +% To get proj.4 string from EPSG, use gdalsrsinfo. Example: +% +% gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" %give ourselves unique file names Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25035) +++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25035) @@ -0,0 +1,55 @@ +from os import remove +import shlex +import subprocess +import tempfile + +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) + + 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) + + 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: + + gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //" + ''' + + # Give ourselves unique file names + file_in = tempfile.NamedTemporaryFile('r+b') + filename_in = file_in.name + file_out = tempfile.NamedTemporaryFile('w+b') + filename_out = file_out.name + + 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)) + subprocess.check_call(args) # Will raise CalledProcessError if return code is not 0 + + A = loadvars(filename_out) + xout = A[0] + xout = xout.reshape(x.shape) + yout = A[1] + yout = yout.reshape(y.shape) + + os.remove(filename_in) + os.remove(filename_out) +#}}} Index: /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py (revision 25035) +++ /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.py (revision 25035) @@ -0,0 +1,38 @@ +import numpy as np + +from mesh3dsurface import mesh3dsurface +from planetradius import planetradius + +def TwoDToThreeD(md, planet): #{{{ + # Reproject model into lat, long if necessary + if md.mesh.proj != epsg2proj(4326): + md.mesh.x, md.mesh.y = gdaltransform(md.mesh.x, md.mesh.y, md.mesh.proj, 'EPSG:4326') + + # Make a 3dsurface mesh out of this + R = planetradius(planet) + + # We assume x and y hold the long, lat values + long = md.mesh.x + lat = md.mesh.y + + # Assume spherical body + 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)) + + elements = md.mesh.elements + vc = md.mesh.vertexconnectivity + vb = md.mesh.vertexonboundary + md.mesh.mesh3dsurface() + md.mesh.lat = lat + md.mesh.long = long + md.mesh.x = x + md.mesh.y = y + md.mesh.z = z + md.mesh.elements = elements + md.mesh.numberofelements = len(elements) + md.mesh.numberofvertices = len(lat) + md.mesh.r = R * np.ones(md.mesh.numberofvertices) + md.mesh.vertexconnectivity = vc + md.mesh.vertexonboundary = vb +#}}} Index: /issm/trunk-jpl/src/m/mesh/meshintersect3d.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshintersect3d.py (revision 25035) +++ /issm/trunk-jpl/src/m/mesh/meshintersect3d.py (revision 25035) @@ -0,0 +1,67 @@ +import math + +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import numpy as np + +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). + ''' + + # Process options + options = pairoptions(*args) + + # Retrieve tolerance + maxtol = options.getfieldvalue('maxtol', 100000) # 100 km + tolincrement = options.getfieldvalue('tolincrement', []) + force = options.getfieldvalue('force', 0) + + # Go through lats, longs and find within tolerance, the index of the corresponding value in lat, long + indices = np.zeros(len(xs)) + + for i in range(len(xs)): + tolerance = 0 + distance = math.sqrt((x - xs[i]) ** 2 + (y - ys[i]) ** 2 + (z - zs[i]) ** 2) + + s = np.where(distance == 0)[0] + if s.size(): + if s.size() > 1: + # We have two vertices that are coincident! Not good + for j in range(len(s)): + plot(x[s[j]], y[s[j]], z[s[j]], c='cyan', s=40) + print('Vertex %i of input mesh coincides with the following ouput mesh vertices ' % i) + print(s) + if force: + indices[i] = s[0] + else: + raise RuntimeError('') + else: + indices[i] = s + else: + # We could not find a 0 distance, find the lowest tolerance that generates a find + count = 1 + while not s.size(): + if count > 1000: + print('could not find a vertex matching vertex %i of input mesh!' % i) + print('Might think anbout changing tolerance increment') + raise RuntimeError('') + tolerance = tolerance + tolincrement + s = np.where(distance < tolerance) + count = count + 1 + if tolerance > maxtol: + print('found matching vertices %i in output mesh for input mesh vertex %i' % (s, i)) + print('however, these vertices are farther than the max tolerance allowed!') + raise RuntimeError('') + + # Recover minimum distance + sf = distance[s] + pos = np.where(sf == sf.min()) + s = s[pos] + indices[i] = s + + if np.where(indices == 0).size(): + raise RuntimeError('issue with transition vector having one empty slot') +#}}}