Other Projections in GSL

If you know the formula for your projection, then it is of course possible to use the data operators within GSL and do the projection yourself. The code provided on this page includes the following projections:

• Orthographic
• Polar stereographic
• Mollweide

function rotate_world(float x0, float y0, float X, float Y,
                              float &X2, float &Y2)
  float cosx, cosy, sinx, siny;
  float i,j,k,ip,jp,kp;

  # The world is rotated so that (x0,y0) goes to the NP
  x0 = x0*pi/180;
  y0 = y0*pi/180;

  cosx = cos(x0); cosy = cos(y0);
  sinx = sin(x0); siny = sin(y0);

  # Convert the data into i,j,k
  i = cos(Y*pi/180)*cos(X*pi/180);
  j = cos(Y*pi/180)*sin(X*pi/180);
  k = sin(Y*pi/180);

  # Rotate i,j,k
  ip = (siny*cosx)*i+(sinx*siny)*j-cosy*k;
  jp = cosx*j-sinx*i;
  kp = (cosx*cosy)*i+(sinx*cosy)*j+siny*k;

  # Go back to A and theta
  Y2 = acos(kp);
  X2 = atan2(jp,ip);

endfunction

function do_orthographic(float x0, float y0, float X, float Y, 
string X2, string Y2)float r;

  rotate_world(x0,y0,X,Y,X,Y);
  r=sin(Y);
  $X2 = if(Y<=pi/2,r*sin(X),undef);
  $Y2 = if(Y<=pi/2,-1*r*cos(X),undef);
endfunction

function do_mercator(float X, float Y, string X2, string Y2)
  $X2=X*pi/180;
  $Y2=log(1 / tan((45-Y/2)*pi/180) );
endfunction

function do_mollweide(float X, float Y, string X2, string Y2)
  $X2=X/90;
  $Y2=cos((90-Y)*pi/180);
endfunction

function do_stereographic(float X, float Y, string X2, string Y2)
  float r;

  r = tan((90-Y)*pi/360);
  $X2 = r*sin(X*pi/180);
  $Y2 = -1*r*cos(X*pi/180);
endfunction

#Example

reset all;
import_worldmap(,,,,,,0,90);

do_orthographic(12,55,Long,Lat,"Long","Lat");
#do_stereographic(Long,Lat,"Long","Lat");
#do_mollweide(Long,Lat,"Long","Lat");

create Viewport page_1.viewport
 ( XuNyRatio = 1
 );
create Domain page_1.viewport.domain;
page_1.viewport.domain.xaxis1.XuNaxle = false;
page_1.viewport.domain.yaxis1.XuNaxle = false;
create Graph page_1.viewport.domain.graph1
 ( XuNgraphType = "line",
   XuNxData = "Long",
   XuNyData = "Lat"
 );