Convert Latitude/Longitude values to x,y on a map

albers equal-area projection

http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html

Sorry, the link for the National Atlas should have been:
http://nationalatlas.gov/natlas/NatlasStart.asp

I've seen the Mathworld site before, but I honestly don't know what to do with this stuff. I'm not a mathematician, I'm a Web developer and artist - so the trigonometry is beyond my capacity. I simply want to find a formula that let's me convert lat/long to x,y on a map of the United States. The size of the map is fixed (600x400) and the lat/long at the corners of the map are known. Given that info I should be able to devise a solution to this problem, it's just a question of coordiate translation and adjusting for curvature of the Earth, right?

here's Visual basic procedure to do it....

you need to specify the  maps origin latitiude and longitude (in radians)  

i got the standard parallels from the map example you have on the site...


hth


Private Sub CalcXY(Lat As Double, Lon As Double, x As Double, y As Double)
' lat you latitude you want
' lng the longitude you want
'
'note pass the lat & lon as radians
'
' radians = degrees * PI / 180
'

Dim mpi as double ' sorry can't remember how to specify pi in vb...
mpi=3.14519


' maplat map origin latitude
' maplon map origin longitude

Dim MapLon As Double
Dim Maplat As Double
' you said you know these ...
maplon = ????? * (MPI / 180)
maplat = ????? * (MPI / 180)


'sp1 standard parallel 1  in your map 28 30' N
'sp2 standard parallel 2  in your map 45  30' N
Dim sp1 As Double
Dim sp2 As Double

' set standard parallels converting degrees to radians
sp1 = 28.5 * (MPI / 180)
sp2 = 45.5 * (MPI / 180)


' Internal calculations
Dim N As Double
Dim C As Double
Dim p0 As Double
Dim O As Double
Dim p As Double

N = 0.5 * (Sin(sp1) + Sin(sp2))
C = (Cos(sp1) * Cos(sp1)) + 2 * N * Sin(sp1)
p0 = Sqr(C - 2 * N * Sin(Maplat)) / N
O = N * (Lon - MapLon)
p = Sqr(C - 2 * N * Sin(Lat))


x = p * Sin(Lat)
y = p0 - p * Cos(Lat)

End Sub

use this for PI

pi=4*atn(1)

Hey Lowfatspread,
Your answer is a bit over my head, but it's clearly thought out and developed. I see the core formula in there, which I'll work on in converting to PHP.
Much thanks.
Jennifer.

Hey everyone,
I've finally got this working using a very simple set of formulas (in PHP). The biggest battle is getting an accurate Mercator map and knowing your map's endpoints. Here's the map solution we used (you can save the output as an Adobe Illustrator file):
http://www.aquarius.geomar.de/omc/make_map.html

I created a map using the above tool with the following endpoints: N: 53, E: -66.54, S: 23, W: -125.49
This creates a nice Continental US & partial Canadian map. Remember to check the "download" box if you want to save the output as a Illustrator file!

I imported the map into Macromedia Flash MX. The map is 600 x 396 (pixels). Here's the working model on my Web site:
http://www.graphiteproductions.com/mapping/flash/map.php

Flash reads a flat text file containing the following result sets (as example):
totDots=10
&nam1=SF State
&loc1=San Francisco, CA
&posx1=33.18
&posy1=226.97
&link1=http://www.sfsu.edu/

&nam2=Wayne State
&loc2=Detroit, MI
&posx2=432.46
&posy2=162.752
&link2=http://www.wayne.edu

&nam3=University of Texas
&loc3=Austin, TX
&posx3=282.34
&posy3=314.289
&link3=http://www.utexas.edu/
...

You get the idea. You could use PHP to obtain this from a database dynamically rather than use a flat text file, if you choose.
Generating the x,y values from long/lat is, of course, the real trick.

Here's how I derived the y coord:
function findYcoord($myLat) {
      $mapHeight = 396;
      $rfactor = 290;
      $radBtm = deg2rad(23);
      $radPixel = deg2rad($myLat);
      $sinRadBtm = sin($radBtm);
      $sinRadPixel = sin($radPixel);
      $convHtBtm = $rfactor * log((1 + $sinRadPixel)/(1 - $sinRadPixel));
      $convHtPixel = $rfactor * log((1 + $sinRadBtm)/(1 - $sinRadBtm));
      $myTotHt = abs($convHtPixel - $convHtBtm);
      $myYcoord = round($mapHeight - $myTotHt, 3);
      return $myYcoord;
}

NOTES:
You simply pass in the latitude value, which is later converted to a radian.
$mapHeight is the height, in pixels, of my map in Flash.
$rfactor is a constant variable. It corresponds to your map scale. You might have to tweak this value up or down a bit (but it worked for me at 290).
$radBtm = deg2rad(23); the value "23" is the latitude at the bottom of my map.

Getting the x value is much, much easier:
$x =round(((125.49 - $lon)*10.178), 2);

NOTES:
125.59 is the longitude at the top left corner of my map.
$lon is the longitude for the city I want to plot
10.178 is a constant that is based on the width of the map, it's calculated as follows:
longitude West - longitude East / width of map
125.49 - 66.54 = 58.95
600 / 58.95 = 10.178
Let me put it another way: every 10.178 pixels on your map equals one degree of longitude.

Here's how I plotted Detroit:
Detroit, MI<br>
Longitude: 83<br>
Latitude: 42.24<br>
<?php
$lon = 83;
$lat = 42.24;
$x =round(((125.49 - $lon)*10.178), 2);
print "x: ".$x."<br>";
print "y: ".findYcoord($lat);      
?>

Works like a champ, I swear by the results and it's far less complicated or confusing than any other solution I've seen thus far.
Thanks again to everyone.
David and Jennifer.

Hi Jennifer,
               Please try this link.i hope this will solve ur problem.
http://www.dmap.co.uk/ll2tm.htm
Ark

Jennifer,

Just curious...how did you calculate the rfactor?

daryn

Good question.
rfactor is an arbitrary value, I just make it up. Trial and error reveals the best results. I know, less than empirical but who said I was perfect?
J.

jloring

Thanks so much for the help. I was researching this exact implementation. If you figure out the rfactor calculation let me know. If I can figure it out I will post.

F

Has anyone figured out the rfactor constant yet?

I've been trying to for a while now - I have implemented this formulae into ASP to project map points from a GPS receiver.

Because I'm using different maps, and switching between them, I the formulae to work out the rfactor, rather than implementing it manually.

Things I've learned about the rfactor:

1.  It is based on the height of your graphic - if your graphic height changes, regardless of the lat/long boundaries, the figure will change.

2.  It is latitude dependant - the larger scale the map is (ie higher difference between top lat and bottom lat) the smaller the number becomes.  Know this as I have a region of the south of the UK and am using an rfactor of 49500 rather than an rfactor of 200 for the whole of the USA!

3.  Possibly latitude location dependant - the figure could be determined by the position of the top and bottom lats based on the vector, ie figure could change based on being closer or further from the equator (0deg lat, no distortion)

I welcome your comments please - jloring, and anyone else who has experience with this, if you are reading this and have a conclusion to this story, please let me know.

Thanks in advance,

Glenn

THIS IS WHAT WORKED FOR ME. I IGNORED THE R FACTOR, AND MADE THE MAP VS. PIXEL CONVERSION DYNAMIC. IT IS COMMENTED FOR YOUR CONVINIENCE. HOPE THIS HELPS YOU ALL OUT. YOUR COMMENTS HELPED ME:
BY THE WAY, THIS IS WRITTEN IN COLDFUSION, BUT IS SIMILAR TO MOST OTHER LANGUAGES.
<!---
MAPVIEW.CFM (c)2004 Trios Associates Inc. Authored by Frank Triantos 6-17-04
this code takes in the parameters of a map in pixels and the latitude and longitude boundaries of said map
allowing a query to be generated that pulls all instances inside said boundaries
and then displays the hot linked position (sing a red-dot image) on top of the map image
which can be clicked to generate another page specific to that item...
----->
<cfscript>
      //enter the height and width of the map image in pixels
      pixel_height = 700;
      pixel_width = 1200;
      // this function converts a lat or long integer to seconds
      function toSeconds (input_a) {
            variables.input_asec = right(input_a, 2);
            variables.input_amin =  mid(input_a, (len(input_a) - 3) , 2) * 60;
            variables.input_adeg =  mid(input_a, (len(input_a) - 6) , 3) * 3600;
            output_a = (input_asec + input_amin + input_adeg);
            return output_a;
      }
      //enter the higher and lower absolute value of longitude and longitude
      long_small = 0590000;
      long_big = 1320000;
      lat_small = 0210000;
      lat_big = 0510000;
      //these values are then converted to seconds for use in setting the map
      //...vertical scale info determined (seconds per pixel)
      long_small_s = toSeconds(long_small);
      long_big_s = toSeconds(long_big);
      secs_width = (long_big_s - long_small_s);
      horizontal_scale = secs_width;
      //...horizontal scale info determined (sedonds per pixel)
      lat_small_s = toSeconds(lat_small);
      lat_big_s = toSeconds(lat_big);
      secs_height = (lat_big_s - lat_small_s);
      vertical_scale = (secs_height / pixel_height);
</cfscript>
<html>
<head>
      <title>Untitled</title>
</head>
<script language="javascript">
      //var precache = new Image(20, 20)
      //precache.src = "dot.jpg"
      function GoThere(loc){
            top.left.location.href="../parent.cfm?ConfigID="+loc;
      }
</script>
<body style="margin:0px;">
<!--------------------------------Query specific to database-------------------------------------------------->
<cfquery name="LetsEditor" datasource="#variables.ds#">
      SELECT intLatitude, intLongitude, configID
      FROM tblConfig WHERE configID IN (2141,4394,19193);
</cfquery>
<!---------------looping through the results of the query----------------------------------------------------->
<cfoutput query="LetsEditor">
<!----scripting section to generate the position of the item in relation to the map parameters---------------->
      <cfscript>
            //verifies that coordinates are in integer format and seperates (from right) seconds, minutes and degrees
            //...converting each to seconds
            //...beginning with x
            
            dd_x = int(intLongitude);
            dd_y = int(intLatitude);
            variables.dd_xsec = right(intLongitude, 2);
            variables.dd_xmin =  mid(intLongitude, (len(intLongitude) - 3) , 2) * 60;
            if (len(intLongitude) EQ 7){
                  variables.dd_xdeg =  mid(intLongitude, (len(intLongitude) - 6) , 3) * 3600;
            } else {
                  variables.dd_xdeg =  mid(intLongitude, (len(intLongitude) - 5) , 2) * 3600;
            }
            //...then for the y
            variables.dd_ysec = right(intLatitude, 2);
            variables.dd_ymin =  mid(intLatitude, (len(intLatitude) - 3) , 2) * 60;
            if (len(intLatitude) EQ 7){
                  variables.dd_ydeg =  mid(intLatitude, (len(intLatitude) - 6) , 3) * 3600;
            } else {
                  variables.dd_ydeg =  mid(intLatitude, (len(intLatitude) - 5) , 2) * 3600;
            }
            //then totaling the converted segments into one whole coordinate (in seconds)
            dd_x = (dd_xsec + dd_xmin + dd_xdeg);
            dd_y = (dd_ysec + dd_ymin + dd_ydeg);
            //pixel position on the map is
            //...for x the left-most (in seconds) coordinate minus the coordinate in question
            //...equals amount from left side multiplied by the quantity
            //...of total mapwidth (pixels) over total of seconds displayed on map
            //...(giving us pixels per second)  
            variables.rel_x = (((long_big_s - dd_x) * pixel_width) / secs_width);
            //for y we have to subtract the conversion from the higher converted value
            //...(which is the lower height on the map because these are negative)
            //...multiply the seconds per pixel conversion
            //...and subtract from the height of the picture in pixels
            variables.rel_y = (((lat_big_s - dd_y) * pixel_height) / secs_height);
      </cfscript>
      <!---------positioning the red dot image using the coordinates located above (scripting section)----------->
      <!---------includes a hot link to the page with more info and passes the id of the item via url------------>
      <div ID="dot#letseditor.configid#" style="position:absolute; left:#variables.rel_x#; Top:#variables.rel_y#;">
            <a href="##" onClick="GoThere(#LetsEditor.ConfigID#)">
                  <img src="dot.jpg" alt="#intLatitude#,#intLongitude#" width="10" height="10" border="0">
            </a>
      </div>
</cfoutput>
<!---insertion of the map image-------------------------------------------------------------------------------->
<img src="usa.gif" name="usa">
</body>
</html>
<!-----------------------------------------ENDING MAPVIEW.CFM-------------------------------------------------->

@ ftriantos:

Firstly, thank you very much for sharing your source code.  I can see the relation between difference in latbig and latsmall and pixels - a straight pixels per lat conversion.

Can we see an online demo of your work?  I would be interested in seeing this....  My ASP isnt on an internet connected server yet, as it is still in testing for accuracy - another reason for the rfactor conversion requirement.

My only comment would be how do you compensate for spherical distortion?  There is an old tutorial in PHP somewhere which uses a world map, with a spheriacally distorted map, so straight x/y to lat/lon conversions are okay, but I'm using a vector of the world, and a map that has been spherically distorted.  jlorings formulae works for me - just this dammed rfactor to decipher!!!

Lowfatspreads answer just killed my brain!!  I can see the workings, just cant put it together - lowfatspread, if you are listening, can you be more specific on the values that need to be passed (ie are the standard parallels the top lat and bottom lat?  The origin lat and lon? is that the top left lat and lon? If so, why do you need the standard parallels in there? How is x converted from just an origin lat and lon?)  I'm sure this is probably very annoying for you, but my maths is extremely weak - I do try and understand, but your understanding outweighs mine tenfold....humble, humble...........

I have found some more info, looking into the trig calcs (although I'm definitely no mathematician!)

1.  rfactor for USA is 290 for USA.
2.  I'm using 49500 for UK
3.  Both are in Northern Hemisphere (above equator - 0 deg lat)
4.  rfactor is higher for a more magnified area with lower distance
5.  rfactor is used in only two areas of conversion - on convHtBtm and convHtPixel

I will try and convert your CFM into ASP and give it a go anyways - it's obviously working well for you and good job on the coding - nice to see comments showing whats going on, takes the pain out of deciphering code!!  Also going back to the drawing board on the rfactor - gonna transopse the formulae and see what happens.......  keep you and all else watching this space posted!!

All,

I have been trying to design a Flash Application that geo-references a vector graphic to enable accurate plotting of "sites of interest".  I have a working proof of concept that seems to be extremely accurate.  Sorry that it is so Flash specific but the general concept of how to do it should work for most circumstances.  I used ColdFusion to query the data base and run the calculations.  If anyone is interested, I can post the CF as well.  Here is a general description of what I did:

First you will need to import the map(vector) graphic you will be using.

Make sure that the graphic fully covers the stage.  (You may have to adjust the size of the stage and or the deminsions of the graphic to accomplish this.)

Select the lat/long intersection that will serve as your max values.  (i.e. for my map I selected 180W and 75N.

Adjust your graphic to be at 0,0 X,Y on the Flash stage at that max lat/long intersection.

Now the fun part:

Select the values for the following that you want your plot points to be between:
maxLong = (in my case 180)
minLong = (in my case 0)
maxLat = (in my case 75)
minLat = (in my case 0)

Measure the distance on the X-axis from 0 pixels to the furthest east longitude line you want to use in your plots.  (i.e. I wanted to plot points that fell between 180w and 0 degrees.  The distance in Flash from the 180w and 0, in pixels was 1138.3.) This gives us a width variable:

xwidth = (in my case 1138.3 px)

Now measure the distance between your upper most latitude line (in my case 75N) and the next latitude line you come to below it (in my case 60n).  You will do this for all your latitude lines.  Measure them from the previous to the next and also note the overall distance from the Y=0 point.  So in my case I have 5 blocks of latitude (from 75-60, 60-45, 45-30, 30-15, and 15-0).  The notation should look like this:

latGroupA:  (75N - 60N)
maxGroupLat = (in my case 75)
minGroupLat = (in my case 60)
yDistanceA = (in my case 256.3px)
yheightA = (in my case 256.3px)

latGroupB:  (60N - 45N)
maxGroupLat = (in my case 60)
minGroupLat = (in my case 45)
yDistanceB = (in my case 412.8px)
yheightB = (in my case 156.5px) or (latGroupB:yDistance - latGroupA:yDistance)

latGroupC:  (45N - 30N)
maxGroupLat = (in my case 45)
minGroupLat = (in my case 30)
yDistanceC = (in my case 532.1px)
yheightC = (in my case 119.3px)

latGroupD:  (30N - 15N)
maxGroupLat = (in my case 30)
minGroupLat = (in my case 15)
yDistanceD = (in my case 634px)
yheightD = (in my case 101.9px)

latGroupE:  (15N - 0)
maxGroupLat = (in my case 15)
minGroupLat = (in my case 0)
yDistanceE = (in my case 728.4px)
yheightE = (in my case 94.4px)

Now, you will need to figure out the pixels per degrees of both your X and Y axis.  The X-axis is fairly easy:

xwidth/(maxLong - minLong) = xFactor

in my case:  1138.3/(180 - 0) = 6.3238888 = xFactor

The Y-axis is a bit more challenging.  You have to determine a pixels per degree for each latGroup.  The equation is virtually the same:

yheight/(maxGroupLat - minGroupLat) = yFactor

So in my case:
latGroupA:
256.3/(75 - 60) = 17.086666 = yFactorA

latGroupB:
156.5/(60 - 45) = 10.433333 = yFactorB

latGroupC:
119.3/(45 - 30) = 7.9533333 = yFactorC

latGroupD:
101.9/(30 - 15) = 6.7933333 = yFactorD

latGroupE:
94.4/(15 - 0) = 6.2933333 = yFactorE

Now to find your X,Y coordiinates:
For X:
(maxLong - ABS(longToPlot)) * xFactor = xCoord

So in my case for Chicago O'Hare International Airport(41.9786, -87.9048):
(180 - ABS(-87.9048)) * 6.3238888 = 582.3998

It's a bit more complex for Y:
First look at the latitude you want to plot and figure out which latGroup it falls within.  Use that groups maxGroupLat and yFactor in the following way:
((maxGroupLat - ABS(latToPlot)) * yFactor) + previous:latGroup:yDistance = yCoord

So in my case for Chicago O'Hare International (41.9786, -87.9048):
((45 - ABS(41.9786)) * 7.9533333) + 412.8 = 436.8302

Hope this helps...

Daryn

 

back @ SixGunZ

http://public.triosinc.com/frank/

By using the mercator map that jLorings comment from 2/8 directed me to (and altering it some) I came up with a mapping procedure without the need to consider spherical distortion. That being said, it is basically only good for the continental US, and I would like to look a little more deeply at the problem. It looks like the comment posted by Daryn above might be the solution I eventually will migrate to. I posted a demo site for you to look at, using three cities that were in my database. These are the three I used to realize that it was working.

I will be continuing on this effort and will post the eventual result when I am finished. What you will see at my link is the exact code I published in the original posting above (mapview.cfm)

Frank