Convert grahpic map X/Y pixels from/to Latitude/longitude

[HSPalm]

Hi all.
I am doing a programming project where I have a map which is 879x436 pixels. At all endpoints I have the geographical coordinates in latitude and longitude corresponding to 0x0, 0x436, 879x0 and 879x436 pixels endpoints X/Y in all corners. My biggest problem is how to calculate the lat/long from a given pixel coordinate.

Say I need the lat/long from the position (0,218) which is half way down the Y axis on the left. If I just divide the latitude difference from Y=0 to Y=436 by 2 I should have the answer, but the problem is that lat/long coordinates are given by angle.minute.seconds. So I could end up with something like 52deg.80(!)min.XXseconds which is a faulty minute.

Don't mind the Earth beeing an ellipsoid, the map is about 7x7 km so I have chosen to ignore that.

My second problem is that latitude on both edges at X axis do not correspond, same goes for longtitude. This may be my map, because I see that the same problem to not occur in google maps. But never mind that just yet :)

I am very greatfull for any contributions!

 

[LCKurtz]

So I could end up with something like 52deg.80(!)min.XXseconds which is a faulty minute.

80 min = 1 degree + 20 min. Just carry the one and call it 53 deg 20 min.

 

[HSPalm]

Yes, thanks, it seems obvious now that you point it out. But its not obvious to me that I don't "miss" any minutes/seconds along the calculation by just dividing by two? Do I?

 

[LCKurtz]

Yes, thanks, it seems obvious now that you point it out. But its not obvious to me that I don't "miss" any minutes/seconds along the calculation by just dividing by two? Do I?

I didn't respond to the rest of your question because I don't have time to research it to be sure I give you an accurate answer. Conversion would depend on what projection was used to map the sphere to a flat surface in the first place. One example is the Mercator projection which you can read about here:

http://en.wikipedia.org/wiki/Mercator_projection

Perhaps that will help you.

 

[HSPalm]

I didn't respond to the rest of your question because I don't have time to research it to be sure I give you an accurate answer. Conversion would depend on what projection was used to map the sphere to a flat surface in the first place. One example is the Mercator projection which you can read about here:

http://en.wikipedia.org/wiki/Mercator_projection

Perhaps that will help you.

Well, maybe this is supposed to tell me something: "WGS84 geographically Lat / Lon (DMM)".
I have read about the Mercator projection, and I wish my map was projected like this, because I don't really need to hassle with how they flattened my map, because the map is so small. I want to assume I can have the same longitude in all Y positions given an X value. But anyway, thank you for your contributions so far, it really means a lot to me.

 

[LCKurtz]

Well, maybe this is supposed to tell me something: "WGS84 geographically Lat / Lon (DMM)".
I have read about the Mercator projection, and I wish my map was projected like this, because I don't really need to hassle with how they flattened my map, because the map is so small. I want to assume I can have the same longitude in all Y positions given an X value. But anyway, thank you for your contributions so far, it really means a lot to me.

If the region is small enough, the function you are seeking would be approximately linear in two variables. Instead of using your pixel numbers, let's call the lower left corner (a,b), the lower right corner (c,d), the upper left corner (f,g), and the upper right corner (k,l). You can put the pixel values in later.

Lets call the (latitude,longitude) values of the lower left corner (u,v), of the lower right corner (r,s), the upper left corner (p,q), and the upper right one (m,n). Then you might assume linear formulas:

X = Ax + By + C
Y = Dx + Ey + F

where you put in the pixel values for (x,y) and get (lat,long) for (X,Y). For example, the lower left corner values would give you

u = Aa + Bb + C
v = Da + Eb + F

Using the lower right and upper left corners listed above will give you 4 more equations. You will have 6 equations in the unknowns A,B,C,D,E,F. These can easily be solved, for example, Maple gives:

A = (-g*u+d*u-d*p-r*b+r*g+p*b)/(f*b-c*b+c*g-a*g+a*d-f*d)
B = -(-f*u+a*p+c*u-c*p-a*r+f*r)/(f*b-c*b+c*g-a*g+a*d-f*d)
C = (-p*c*b+d*a*p-d*f*u+r*f*b+g*c*u-g*a*r)/(f*b-c*b+c*g-a*g+a*d-f*d)
D = (-g*v+d*v-d*q-s*b+s*g+q*b)/(f*b-c*b+c*g-a*g+a*d-f*d)
E = -(-f*v+a*q+c*v-c*q-a*s+f*s)/(f*b-c*b+c*g-a*g+a*d-f*d)
F = (-q*c*b+d*a*q-d*f*v+s*f*b+g*c*v-g*a*s)/(f*b-c*b+c*g-a*g+a*d-f*d)

So here's what to try. Express your latitude, longitude values either all in terms of decimal degrees or in terms of minutes. Use them and your pixel values to calculate the right side of the above equations. Once you know A,B,C,D,E,F you have the equations:

X = Ax + By + C
Y = Dx + Ey + F

completely known. To see if this is going to work for you, you can now plug the upper right corner (k,l) into these equations and see how close to the upper right values (m,n) you have for (lat,long) you get. If it's close enough, just program in these last two decimal equations and you are home free.

 

[HSPalm]

upper left upper right
pixels(f,g = 0,0) pixels(k,l = 879,0)
lat/lon(p,q = 59 53.2194 , 10 27.321) lat/lon(m,n = 59 53.5566 , 10 37.2594)


lower left lower right
pixels(a,b = 0,436) pixels(c,d = 879,436)
lat/lon(u,v = 59 50.7486 , 10 27.6582) lat/lon(r,s = 59 51.0846 , 10 37.5852)

If I want to express my lat/lon as degrees only, it should be written as deg,mmssss because the minutes and seconds are the "decimal points" of the degrees, right?
So here it goes

A = (25945,2639-25956,03658-25946.72886+25956,03658)/(-383244)=0,00000382252
B = -(52307,08019-52328.79853)/(-383244)=-0,00005666974
C = (-59,532194*879*436)/(-879*436)=59,532194
D = (4480,589752-4479,11956-4523,871472+4479,11956)/(-383244)=0,00011293515
E = -(879*10,27-879*10,27321)/(-879*436)=-0,00000736238
F = (-10,27321*879*436)/(-879*436)=10,27321


so
m = Ak + Bl + C = 0,00000382252*879 + 0,00011293515 = 0,00347293
n = Dk + El + F = 0,00011293515*879 + 10,27321 = 10,273322

n is of by 0.1, but m is way of. Arggh, this took my at least a couple of hours. Hmm. I see how all of this is right, I don't see what I have done wrong. Guess I should have written down all my steps in my calculation.

But one more thing. Is this method depending on the latitude beeing the same in all x (pixel) positions? Because mine are not, they are slightly different, which I guess is natural. Maybe I should use the mean of p and m for both p and m, the mean of u and r for both u and r and so on to get a rectangular projection of my map.

edit:
I just realized something... lat and long are given by "deg min.ssss", but why do my gps output latitudes like 59 52.786? There are only 60 seconds in one minute, so why do I get a reading of 52 minutes 78 seconds? SO, since my map spans only between minutes, and not degrees (which actually changes once per 60 minutes, not per 100 like the seconds) I do not have to think about converting between 0-100 and 0-60 values! This is great.
So if upper left is (0,0) pixels and 59.532194 latitude
and lower left is (0,436) pixels and 59.507486 latitude
and I want to know the latitude halfway between pixel y=0 and y=436 (which is 218) at x=0
then I know: latitude=(Yposition/436)*(59.532194-59.507486)+59.507486=(218/436)*(59.532194-59.507486)+59.507486=59.51984

Am I right?

 

[LCKurtz]

I will probably have some time today to look at this again. I will plug your numbers into Maple to see if I agree with what you get.

 

[HSPalm]

So that's what maple is, software? Sorry, I had no clue!

 

[LCKurtz]

So that's what maple is, software? Sorry, I had no clue!

Yes, Maple is a powerful software package. I have put your numbers in. If you express the latitude and longitude figures in minutes you get:

(a,b) = (0,436), (u,v) = (3590.7486,627.6582)
(c,d) = (879,436), (r,s) = (3591.0846,637.5852)
(f,g)=(0,0), (p,q)=(3593.2194,627.321)

These give these two equations for X and Y:

X = 0.0003822525597 x - 0.005666972477 y + 3593.219400
Y = 0.01129351536 x + 0.0007733944954 y + 627.3210000

When put the fourth point (k,l) = (879,0) in for x and y to calculate the fourth point you get

(m,n) = (3593.555400,637.2480000). These are expressed in minutes, and converting them to degrees and minutes gives

m = 59 deg 53.555400 min
n = 10 deg 37.2480000 min.

Pretty close, actually, so those two equations for X and Y may be what you need.

 

[HSPalm]

Yes, Maple is a powerful software package. I have put your numbers in. If you express the latitude and longitude figures in minutes you get:

(a,b) = (0,436), (u,v) = (3590.7486,627.6582)
(c,d) = (879,436), (r,s) = (3591.0846,637.5852)
(f,g)=(0,0), (p,q)=(3593.2194,627.321)

These give these two equations for X and Y:

X = 0.0003822525597 x - 0.005666972477 y + 3593.219400
Y = 0.01129351536 x + 0.0007733944954 y + 627.3210000

When put the fourth point (k,l) = (879,0) in for x and y to calculate the fourth point you get

(m,n) = (3593.555400,637.2480000). These are expressed in minutes, and converting them to degrees and minutes gives

m = 59 deg 53.555400 min
n = 10 deg 37.2480000 min.

Pretty close, actually, so those two equations for X and Y may be what you need.

Thank you so much! The accuracy is more than enough for me. But how did you convert between minutes and degrees? Especially between minutes and degrees,minutes?

I don't want to be a jerk and say it wasn't worth anything, but do you agree that the whole conversion is not necessary in my specific situation? Because the span of my map is so small it doesn't exceed a whole hour, and the seconds in latitude and longitude are not really seconds, but just decimal points?

 

[LCKurtz]

Thank you so much! The accuracy is more than enough for me. But how did you convert between minutes and degrees? Especially between minutes and degrees,minutes?

I don't want to be a jerk and say it wasn't worth anything, but do you agree that the whole conversion is not necessary in my specific situation? Because the span of my map is so small it doesn't exceed a whole hour, and the seconds in latitude and longitude are not really seconds, but just decimal points?

Take 3593.555400 minutes. There are 60 minutes in a degree, so divide that by 60 to see how many degrees there are. This gives 59.89259 degrees. So you have 59 degrees and .89259 part of a degree. Since 1 degree is 60 minutes, to see how many minutes .89259 of a degree is just multiply .89259 * 60 = 53.5554 minutes, which is the form in which your numbers were given. That's about 53 and a half minutes and you could similarly convert the .5554 of a minute to seconds if you wanted.

 

[Sone]

Dear LCKurtz,

Your solution is really perfect for the calculation of Lat/Long from X/Y.
But could you help me to calculate in the other way round? I mean calculate X/Y from Lat/Long while keeping the same input as above (a, b) (c, d) ... (u, v) (r, s )... and same accuracy as above is also fine.

Thank you for your help

 

[LCKurtz]

Dear LCKurtz,

Your solution is really perfect for the calculation of Lat/Long from X/Y.
But could you help me to calculate in the other way round? I mean calculate X/Y from Lat/Long while keeping the same input as above (a, b) (c, d) ... (u, v) (r, s )... and same accuracy as above is also fine.

Thank you for your help

I will give you some numbers this weekend. Got too much company visiting right now to spend any time...

 

[Sone]

Many Thanks and take your time.

 

[LCKurtz]

Dear LCKurtz,

Your solution is really perfect for the calculation of Lat/Long from X/Y.
But could you help me to calculate in the other way round? I mean calculate X/Y from Lat/Long while keeping the same input as above (a, b) (c, d) ... (u, v) (r, s )... and same accuracy as above is also fine.

Thank you for your help

I will give you some numbers this weekend. Got too much company visiting right now to spend any time...

OK, I have some equations for you:

xpixels = 12.02871768*LAT + 88.13925162 LONG - 98513.42519
ypixels = 5.945229958*LONG + 6.274183271*10⁵ - 175.6496958*LAT

where LAT and LONG must be expressed in minutes and you will want to round off the xpixels and ypixels answers to the nearest whole number.

 

[sultanp]

Hi - Thanks very much for the previous work on this question.

I am wondering if it is possible to use all four corners to generate the
linear equations using MAPLE? The above example uses only three,
producing 6 equations. Is there a tighter solution that would use four
corners and eight equations?

Thanks! - Peter

 

[LCKurtz]

Hi - Thanks very much for the previous work on this question.

I am wondering if it is possible to use all four corners to generate the
linear equations using MAPLE? The above example uses only three,
producing 6 equations. Is there a tighter solution that would use four
corners and eight equations?

Thanks! - Peter

No. This solution uses the approximation that the spherical Earth is essentially flat over a small region (we are ignoring topography here) and is approximated well by a region of a plane. It only takes three noncollinear points to determine the equation of a plane. If you try to add a fourth you will either get the same values if the fourth point is exactly coplanar with the others or an inconsistent system if the four points aren't actually quite coplanar.

 

[sultanp]

Thanks.

Wondering if there is a '4-corner' solution for larger areas (several km to
a side) when the flat Earth approximation doesn't hold.

 

[Stephen Tashi]

Well, maybe this is supposed to tell me something: "WGS84 geographically Lat / Lon (DMM)".

The location of a point on the Earth in UTM (Universal Transverse Mecator) coordinates depends on what shape is assumed for the earth. Also, you can't exactly determine a UTM coordinate from Lat-Lon alone since the elevation of the point matters. The projection from a point in space perpendicular to the surface of the Earth is not a projection onto a perfect sphere. There are various models for the non-spherical shape of the earth. These are called "spheroids" and WGS84 is a particular spheroid, commonly used on military maps.

I don't know whether your program is a "serious" or recreational endeavor. If you are dealing with a "serious" application, you'll have to study all the non-intutive intricacies of map coordinate systems. One irritating feature of UTM coordinate maps is that maps of adjacent areas on the Earth sometimes have a "seam" in them when you put them side by side. This is because UTM coordinates use projection on a plane, but the plane used depends on the general area of the Earth you are working in. If you get two adjacent areas that use different planes, the maps don't match at a common boundary.