Contents‎ > ‎

Azimuthal Projection


as of 140619: the program was updated for .Net framework 4.5 with some minor updates

Azimuthal equidistant projection is a way to display the surface of the Earth (or other spherical bodies). A projection of a sphere will always be distorted because the curved surface is mapped to a flat plane. In order to display a map as accurate as possible, it is necessary to choose which information to discard and which to keep. This projection preserves the distance and the direction from the center of the map; however, area is not preserved. For example, in the image above, the areas of Africa and Australia appear the same when in fact this is not the case. The degree of distortion depends on the position in the projection. The distortion at the center is less than that in the peripheral regions.


I made the program to demonstrate this kind of map.
To see the javascript version, go to http://sites.google.com/site/somedayuniverse/contents/azimap_js.

In addition, I have provided an executable version on your PC. Click on this image to download the compressed file.
*Note: this executable requires .Net framework. Please keep its version latest.


The following sections explain the method to create this application. I use a plate carree image as the original format.



To generate this projection, the computer must map each pixel to the correct position. The original data will be taken from the projection called "plate carree." The position of a pixel in plate carree represents its latitude and longitude. When the computer calculates the longitude and latitude of the desired pixel in the projection, it can then use them to 'look-up' the corresponding pixel on the original data. Thus, it is necessary to find the algorithm to convert the position on azimuthal equidistant projection to the location on plate carree.



First we need to list the parameters which are known and which ones need to be calculated.

Variables which are already known:
    the coordinates of the center of the projection on plate carree
    the angular distance between the center and the north pole (derived from "the actual distance / the radius of Earth")
    the angular distance between the center and the target point (derived from "PI/2 - the latitude of Center of Map"), and
    the direction from the north pole to the target point

Items to find:
    the location of the target point on plate carree

This is quite a complex problem to solve. Some simplification is needed to make it easier to visualize.



In this case, the relationship between the points can be simplified into a tetrahedron. In the above image, N is the north pole; A is the center of the Map; B is the target point; and C is the center of the Earth. The condition is so specific that three of its lines are the same as the radius of the Earth.

We need to solve for angle NCB, (which is "PI/2 - B's latitude"), and the angle between plane NCA and plane NCB, (which is the difference in longitude of A and B).


First we will find angle NCB:

Lines NC and BC are radii, so it is easy to solve for angle NCB (by using the cosine rule) if you know the length of NB. Let's say there are additional points, N', B', P, and Q, in the above image. N'N and BB' are parallel to AC. BP, B'Q, N'P, and NQ are perpendicular to AC. In addition, angle N'PB and angle NQB' are the same as the angle between the north pole and the target point, (this angle is represented as theta or θ in the following equations). Therefore, there are the following relationships:



            *R is the radius of Earth. Dist is the distance from A to B. LatA is the latitude of A.
            **this method is based on the law of cosines.

You can find the length of NB using the Pythagorean theorem, and angle NCB.



You finally know the latitude of the target point when you subtract the angle NCB from PI/2. When the value is negative, it means it is the south latitude.



The angle between plane NAC and plane NBC, (which equals angle APB' and angle A'QB in the above) can be found by calculating its cosine. So the three lines, AP, B'P, and AB', are required. AP and B'P are easy to find and AB' can be derived from AB and AA' using Pythagorean theorem.



You then know the longitude of the target point when you add angle APB' to the longitude of the center of the map.


Programming
The following is a part of some sample code. It is written in C++/CLI.

Main variables in the process:
*Some of the definitions depend on your IDE and other environments.
    array<Byte,2>^ sdt; //the source image data
    int sw; //width of the source image
    int sh; //height of the source image
    int sx; //horizontal index (0 <= sx < sw)
    int sy; //vertical index (0 <= sy < sh)

   
array<Byte,2>^ ddt; //the data of the desired map
    int dw; //width of the map
    int dh; //height of the map
    int dx; //horizontal index (0 <= dx < dw)
    int dy; //vertical index (0 <= dy < dh)

    double ltA; //PI/2 - Latitude of A (Center of Map)
    double lgA; //Longitude of A (Center of Map)
    double ltB; //PI/2 - Latitude of B (Target Point)
    double lgB; //Longitude of B (Target Point)

    double leng; //Length from Center of Map to Target Point on the desired map
    double rad; //Angular distance from Center of Map to Target Point (rad = leng / radius of the map * PI)

  Examples of the coordinates
(ltA and lgA)
    If you want the map whose center is Rome, Italy (N41.9 E12.5), the latitude (ltA) and the longitude (lgA) are...
        double ltA = (90 - 41.9) / 180 * PI;
        double lgA = (12.5) / 180 * PI;
    In case of Rio De Janeiro, Brazil (S22.9 W43.2),
        double ltA = (90 - (-22.9)) / 180 * PI;
        double lgA = (- 43.2) / 180 * PI;

MACRO
When I coded this method in C++, I wrote it as a macro to optimise program speed as a macro doesn't have the overhead of a function call. I wrote the macro as the following.

Before calling this, you need calculate these parameters.
    bool dirE: true when the target is in east of the map center
    double sltA: sine of ltA,   e.g. north pole: sltA=0, equator: sltA=1, south pole: sltA=0
    double cltA: cosine of ltA,   e.g. north pole: sltA=1, equator: sltA=0, south pole: sltA=-1
    double cdir: cosine of (the direction from the north pole to B)
                      The azimuth can be derived from acos((dh / 2 - dy) / leng), i.e. cdir = (dh / 2 - dy) / leng
    double sdr: sine of (the angular distance to Target Point), i.e. sdr = sin(rad)
    double cdr: cosine of (the angular distance to Target Point), i.e. cdr = cos(rad)

After the macro, the following variables have been updated:
    double ltB: PI/2 - the latitude of B
    double dlgB: the difference between the longitude A and B

By using dlgB, you also know lgB (= lgA + dlgB)

====================================================================
#define PI 3.14159265358979323846

#define azi_calc(dirE,sltA,cltA,cdir,sdr,cdr,ltB,dlgB) \
    do {\
        /* ltB */\
        double cltB = sltA*sdr*cdir+cltA*cdr;\
        if(cltB > 1)\
            cltB = 1;\
        else if(cltB < -1)\
            cltB = -1;\
        ltB = acos(cltB);\
        /* dlgB */\
        if(cltA == 1){/* when A is on the north pole */\
            if(dirE)\
                dlgB = PI - acos(cdir);\
            else\
                dlgB = acos(cdir) - PI;\
        }else if(cltA == -1){/* when A is on the south pole */\
            if(dirE)\
                dlgB = acos(cdir);\
            else\
                dlgB = -1*acos(cdir);\
        }else if(cdir == 1){ /* when B is in the true north */\
            if(sltA*cdr-cltA*sdr >= 0)\
                dlgB = 0;\
            else\
                dlgB = PI;\
        }else if(cdir == -1){ /* when B is in the true south */\
            if(sltA*cdr+cltA*sdr >= 0)\
                dlgB = 0;\
            else\
                dlgB = PI;\
        }else{ /* others */\
            double sltB = sqrt(1-cltB*cltB);\
            double buf = (cdr-cltA*cltB)/(sltA*sltB);\
            if(buf > 1)\
                buf = 1;\
            else if(buf < -1)\
                buf = -1;\
            if(dirE)\
                dlgB = acos(buf);\
            else\
                dlgB = (-1)*acos(buf);\
        }\
    }while(0)
====================================================================

IMPLEMENTATION

Precalculation of cdir,sdr,cdr
Each pixel of the desired map has specific values for cdir, sdr, and cdr. If you calculate these data in advance, you can reduce time to calculate in the main process.
====================================================================
    array<double,2>^ tbl_cdir;
    array<double,2>^ tbl_sdr;
    array<double,2>^ tbl_cdr;

    tbl_cdir = gcnew array<double,2>(dw,dh);
    tbl_sdr = gcnew array<double,2>(dw,dh);
    tbl_cdr = gcnew array<double,2>(dw,dh);

    int dx = 0;
    while(dx < dw){
        int dy = 0;
        while(dy < dh){
            double leng = sqrt(pow(dw/2-dx,2)+pow(dh/2-dy,2)); //the length from Center of Map to Target Point on the map
            double rad = leng/(dw/2)*PI; //the angular distance from Center of Map to Target Point

            if(leng > dw/2){ //out of the circle of the map
               
tbl_cdir[dx,dy] = -1;
               
tbl_sdr[dx,dy] = -1;
               
tbl_cdr[dx,dy] = -1;
            }else{
               
tbl_sdr[dx,dy] = sin(rad);
               
tbl_cdr[dx,dy] = cos(rad);
                if(leng == 0)
                   
tbl_cdir[dx,dy] = 1;
                else
                   
tbl_cdir[dx,dy] = (dh/2-dy)/leng;
             }
             dy++;
         }
         dx++;

    }
====================================================================

Main Calculation
In this process, you can just pick the data from the tables of cdir, sdr, and cdr; calculate the coordinates of Target Position; and take the pixel from the source to the desired map.
====================================================================
    double sltA = sin(ltA);
    double cltA = cos(ltA);

    bool dirE;
    double cdir;
    double sdr;
    dboule cdr;
    double ltB;
    double dlgB;
    double lgB;

    int dx = 0;
    while(dx < dw){
        int dy = 0;
        while(dy < dh){
            if(
tbl_sdr[dx,dy] == -1){
                ddt[(int)(dx*3+dy*3*dw)] = (Byte)0; //R
                ddt[(int)(dx*3+dy*3*dw)+1] = (Byte)0; //G
                ddt[(int)(dx*3+dy*3*dw)+2] = (Byte)0; //B
             }else{
                 if(dx >= (int)(dw/2))
                     dirE = true;
                 else
                     dirE = false;
                 cdir =
tbl_cdir[dx,dy];
                 sdr =
tbl_sdr[dx,dy];
                 cdr =
tbl_cdr[dx,dy];
                 azi_calc(dirE,sltA,cltA,cdir,sdr,cdr,ltB,dlgB);    //execution of the macro
                 lgB = dlgB + lgA;
                 if(lgB >= PI)
                     lgB = lgB - 2*PI;
                 else if(lgB < -1*PI)
                     lgB = lgB + 2*PI;
                 sx = (int)((lgB+PI)/(2*PI)*(sw-1));
                 sy = (int)(lgB/PI*(sh-1));

                 ddt[(int)(dx*3+dy*3*dw)] = sdt[(int)(sx*3+sy*3*sw)]; //R
                 ddt[(int)(dx*3+dy*3*dw)+1] = sdt[(int)(sx*3+sy*3*sw)+1]; //G
                 ddt[(int)(dx*3+dy*3*dw)+2] = sdt[(int)(sx*3+sy*3*sw)+2]; //B
             }
             dy++;
         }
         dx++;
     }

====================================================================
   
I also used some other techniques and functions such as Marshal::Copy from .Net framework to reduce the calculation time. They might have reduced more time of calculation.
Anyway, if you do not understand my poor writing, please let me know.


*Note A friend from Australia edited this page and helped me correcting my wrong sentences on this page. Thanks a lot!


written on 12/03/10
updated on 12/07/24
updated on 12/08/03
updated on 14/06/19