Changeset 20111

Timestamp:

02/10/16 12:35:12 (9 years ago)

Author:

Eric.Larour

Message:

CHG: new improved patchglobe that relies on gdaltransform instead of ll2xy. Also fixed large issue with connectivity
not being correctly handled when augmenting the mesh with patches.

File:

• issm/trunk-jpl/src/m/mesh/patchglobe.m (modified) (5 diffs)

issm/trunk-jpl/src/m/mesh/patchglobe.m

@@ -5 +5 @@
 
         %recover basic options:
-        hem=getfieldvalue(options,'hem');
         bandwidth=getfieldvalue(options,'bandwidth',100000);
 
         %give ourselves a unique temporary directory: 
         temproot=tempname; mkdir(temproot);
+
+        %align external segments on 2d model: 
+        mh2d.segments=alignsegments(mh2d.segments);
 
         %figure out domain outline: 
@@ -18 +20 @@
         expcontract([temproot '/PatchBroad.exp'],[temproot '/PatchBroad.exp'],-bandwidth);
 
-        %now flag vertices that are within the contour:
-        [x,y]=ll2xy(mh.lat,mh.long,hem);
-        if hem==1,
-                flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat>=0;
-        else
-                flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat<=0;
-        end
+        %now flag vertices (from mh2d's broad contour that are on the global mesh: do this in the local 2d mesh reference system. 
+        [x,y]=gdaltransform(mh.long,mh.lat,'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
+        flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1);
 
         %expand flags to any element that touches the contour: 
@@ -37 +35 @@
 
         %x,y for segments: 
-        [xsegs,ysegs]=ll2xy(mh.lat(mh.segments(:,1)),mh.long(mh.segments(:,1)),hem);
+        [xsegs,ysegs]=gdaltransform(mh.long(mh.segments(:,1)),mh.lat(mh.segments(:,1)),'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
 
-        %create contour out of these segments:
+        %create lat,long contour out of these segments:
         meshtodomain(mh,[temproot '/PatchEnveloppe.exp'],'latlong','on');
-        expll2xy([temproot '/PatchEnveloppe.exp'],hem);
+                
+        %get these lat,long transformed to local mesh referencial:
+        env=expread([temproot '/PatchEnveloppe.exp']); 
+        [env.x,env.y]=gdaltransform(env.x,env.y,'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
 
         %now, create domain outine from broad enveloppe and initial mesh 
-        env=expread([temproot '/PatchEnveloppe.exp']); 
         dom=expread([temproot '/Patch.exp']);
         
@@ -61 +61 @@
         expwrite(domain,[temproot '/PatchBand.exp']);
 
+        %plot mesh: 
+        if  getfieldvalue(options,'plot',0), expdisp([temproot '/PatchBand.exp']); end
+
         %mesh: 
         mdb=bamg(model(),'domain',[temproot '/PatchBand.exp'],'MaxCornerAngle',1e-15,'gradation',10000); 
@@ -71 +74 @@
 
         %augment patch with band
-        [mhb.lat,mhb.long]=xy2ll(mhb.x,mhb.y,hem);
+        [mhb.long,mhb.lat]=gdaltransform(mhb.x,mhb.y,['EPSG:' num2str(mh2d.epsg)],'EPSG:4326');
         mh2db=augment2dmesh(mh2d,mhb);