felixrulz

07-18-2012, 10:00 PM

Hey, this is my first time here... was looking for some help from some cartography savy people

I have created an almost spherical disk by joining lots of triangular/rectangular polygons together and have a equidistant image of the earth loaded as a texture map. What I am having trouble doing is converting/mapping the azimuthal equidistant points on the disk to a latitude and longitude so I can determing the texture coordinate. Here is a snippet of code where the disk is created:

void CreateDisk (double R, double x0, double y0, double z0)

{

const double rad_inc = R/(disk_r_divs);

const double theta_inc = 2*PI/(disk_az_divs);

double theta, radius, b2[4], a2[4];

double latitude, longitude;

int n=0,i; // Current vertex

for(radius=0; radius < (R+rad_inc/2); radius+=rad_inc)

{

b2[0] = radius;

b2[1] = radius + rad_inc;

b2[2] = radius;

b2[3] = radius + rad_inc;

for(theta=0; theta <= PI2; theta+=theta_inc)

{

a2[0] = theta;

a2[1] = theta;

a2[2] = theta + theta_inc;

a2[3] = theta + theta_inc;

for (i=0; i<4; i++)

{

// disk coordinates

p_d[n].x = b2[i] * cos(a2[i]) + x0;

p_d[n].y = b2[i] * sin(a2[i]) + y0;

p_d[n].z = z0;

// TEXTURE COORDINATES CALCULATED HERE

n++;

}

}

}

}

The use of b2[] and a2[] are to calculate the other points of the rectangle. And are not really relavent to the problem I am havving (I think) so this code can be interpreted like this for simplicity:

void CreateDisk (double R, double x0, double y0, double z0)

{

const double rad_inc = R/(disk_r_divs);

const double theta_inc = 2*PI/(disk_az_divs);

double theta, radius;

double latitude, longitude; //used to calculate texture coordinates

int n=0 // Current vertex

for(radius=0; radius < R; radius+=rad_inc)

{

for(theta=0; theta <= PI2; theta+=theta_inc)

{

// disk coordinates

p_d[n].x = b2[i] * cos(a2[i]) + x0;

p_d[n].y = b2[i] * sin(a2[i]) + y0;

p_d[n].z = z0;

// TEXTURE COORDINATES CALCULATED HERE

n++;

}

}

}

This code starts with a radius of 0 and then will sweep through theta (the azimuth of the map). Repeating with larger radii untill reaching R. Here is an example of the resulting disk:

http://dl.dropbox.com/u/49415695/Disk.png

The problem I am havving is converting the cartesian (or polar if this is easier) coordinates to a latitude and longitude for the equirectangular map (This needs to be inserted where indicated in the code sample). I have managed managed a polar aspect by simply using a ratio of the radius and azimuth. (Note that the +0.5 and /PI2 and /R are to convert/normalise the latitude and longitude into a 0:1 range which is required by OpenGL.

p_d[n].u = 0.5 + a2[i]/PI2; // + rotation

p_d[n].v = (b2[i])/R;

This is the resulting image:

http://dl.dropbox.com/u/49415695/Polar%20aspect.png

For the oblique projection I am struggeling. The latest equations I have been trying to use where taken from Map Projections – A Working Manual (pg 196) http://books.google.com.au/books/reader?id=nPdOAAAAMAAJ&printsec=frontcover&output=reader&pg=GBS.PR1,

the equations are reproduced here at Wolfram http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html.

I implemented them in code as follows. b2[] can simply be considered a radius and a2[] an azimuth, p_d[n].y/x are the x/y coordinates:

if (radius == 0)

{

latitude = lat0;

longitude = long0 + theta;

}

else

{

latitude = asin( cos(b2[i])*sin(lat0) + (p_d[n].y*sin(b2[i])*cos(lat0)/b2[i]) );

if (lat0 == PI_2)

longitude = long0 + atan( -p_d[n].x/p_d[n].y );

else if (lat0 == -PI_2)

longitude = long0 + atan( p_d[n].x/p_d[n].y );

else

longitude = long0 + atan( p_d[n].x*sin(b2[i])/( b2[i]*cos(lat0)*cos(b2[i])-p_d[n].y*sin(lat0)*sin(b2[i]) ) );

}

// convert lat/long to texture coordinates

p_d[n].u = longitude/PI2;

p_d[n].v = 0.5 + latitude/PI;

The results seem somewhat correct but I don't even get the entire azimuthal projection (for example with australia at the centre you can't see much else). I have been trying to fix this for near a week and thought I should ask for help.

Thanks.

I have created an almost spherical disk by joining lots of triangular/rectangular polygons together and have a equidistant image of the earth loaded as a texture map. What I am having trouble doing is converting/mapping the azimuthal equidistant points on the disk to a latitude and longitude so I can determing the texture coordinate. Here is a snippet of code where the disk is created:

void CreateDisk (double R, double x0, double y0, double z0)

{

const double rad_inc = R/(disk_r_divs);

const double theta_inc = 2*PI/(disk_az_divs);

double theta, radius, b2[4], a2[4];

double latitude, longitude;

int n=0,i; // Current vertex

for(radius=0; radius < (R+rad_inc/2); radius+=rad_inc)

{

b2[0] = radius;

b2[1] = radius + rad_inc;

b2[2] = radius;

b2[3] = radius + rad_inc;

for(theta=0; theta <= PI2; theta+=theta_inc)

{

a2[0] = theta;

a2[1] = theta;

a2[2] = theta + theta_inc;

a2[3] = theta + theta_inc;

for (i=0; i<4; i++)

{

// disk coordinates

p_d[n].x = b2[i] * cos(a2[i]) + x0;

p_d[n].y = b2[i] * sin(a2[i]) + y0;

p_d[n].z = z0;

// TEXTURE COORDINATES CALCULATED HERE

n++;

}

}

}

}

The use of b2[] and a2[] are to calculate the other points of the rectangle. And are not really relavent to the problem I am havving (I think) so this code can be interpreted like this for simplicity:

void CreateDisk (double R, double x0, double y0, double z0)

{

const double rad_inc = R/(disk_r_divs);

const double theta_inc = 2*PI/(disk_az_divs);

double theta, radius;

double latitude, longitude; //used to calculate texture coordinates

int n=0 // Current vertex

for(radius=0; radius < R; radius+=rad_inc)

{

for(theta=0; theta <= PI2; theta+=theta_inc)

{

// disk coordinates

p_d[n].x = b2[i] * cos(a2[i]) + x0;

p_d[n].y = b2[i] * sin(a2[i]) + y0;

p_d[n].z = z0;

// TEXTURE COORDINATES CALCULATED HERE

n++;

}

}

}

This code starts with a radius of 0 and then will sweep through theta (the azimuth of the map). Repeating with larger radii untill reaching R. Here is an example of the resulting disk:

http://dl.dropbox.com/u/49415695/Disk.png

The problem I am havving is converting the cartesian (or polar if this is easier) coordinates to a latitude and longitude for the equirectangular map (This needs to be inserted where indicated in the code sample). I have managed managed a polar aspect by simply using a ratio of the radius and azimuth. (Note that the +0.5 and /PI2 and /R are to convert/normalise the latitude and longitude into a 0:1 range which is required by OpenGL.

p_d[n].u = 0.5 + a2[i]/PI2; // + rotation

p_d[n].v = (b2[i])/R;

This is the resulting image:

http://dl.dropbox.com/u/49415695/Polar%20aspect.png

For the oblique projection I am struggeling. The latest equations I have been trying to use where taken from Map Projections – A Working Manual (pg 196) http://books.google.com.au/books/reader?id=nPdOAAAAMAAJ&printsec=frontcover&output=reader&pg=GBS.PR1,

the equations are reproduced here at Wolfram http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html.

I implemented them in code as follows. b2[] can simply be considered a radius and a2[] an azimuth, p_d[n].y/x are the x/y coordinates:

if (radius == 0)

{

latitude = lat0;

longitude = long0 + theta;

}

else

{

latitude = asin( cos(b2[i])*sin(lat0) + (p_d[n].y*sin(b2[i])*cos(lat0)/b2[i]) );

if (lat0 == PI_2)

longitude = long0 + atan( -p_d[n].x/p_d[n].y );

else if (lat0 == -PI_2)

longitude = long0 + atan( p_d[n].x/p_d[n].y );

else

longitude = long0 + atan( p_d[n].x*sin(b2[i])/( b2[i]*cos(lat0)*cos(b2[i])-p_d[n].y*sin(lat0)*sin(b2[i]) ) );

}

// convert lat/long to texture coordinates

p_d[n].u = longitude/PI2;

p_d[n].v = 0.5 + latitude/PI;

The results seem somewhat correct but I don't even get the entire azimuthal projection (for example with australia at the centre you can't see much else). I have been trying to fix this for near a week and thought I should ask for help.

Thanks.