Changeset 20523

Timestamp:

04/21/16 14:11:52 (9 years ago)

Author:

Mathieu Morlighem

Message:

CHG: faster conversion and more complete documentation

File:

• issm/trunk-jpl/src/m/coordsystems/gdaltransform.m (modified) (3 diffs)

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

@@ -6 +6 @@
 %
 %   Example: 
-%      [x,y] = gdaltransform(md.mesh.lat,md.mesh.long,'EPSG:3184','EPSG:3411'); 
+%      [x,y] = gdaltransform(md.mesh.long,md.mesh.lat,'EPSG:4326','EPSG:3031')
 %
 %   For reference: 
-%       EPSG: 4326 (lat,long)
-%       EPSG: 3411  (greenland, +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs)
-%       EPSG: 3031 (antarctica, +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs)
+%      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
+%
+%      To get proj.4 string from EPSG, use gdalsrsinfo. Example:    gdalsrsinfo "EPSG:4326" | grep "PROJ.4"
 
         %give ourselves a unique temporary directory: 
@@ -17 +26 @@
 
         fid=fopen([temproot '/.rand1234.txt'],'w');
-        for i=1:length(x),
-                fprintf(fid,'%g %g\n',x(i),y(i));
-        end
+        fprintf(fid,'%8g %8g\n',[x(:) y(:)]');
         fclose(fid);
 
@@ -26 +33 @@
                 error(r);
         end
-        A=textread([temproot '/.rand1235.txt']);
-        xout=A(:,1);
-        yout=A(:,2);
+        A=load([temproot '/.rand1235.txt']);
+        xout=A(:,1); xout=reshape(xout,size(x));
+        yout=A(:,2); yout=reshape(yout,size(y));