5

This had been asked previously on Stack Overflow.

I've tried different methods/formulas but results are not correct. 40.764296,-73.97302 should result into 930,3590 pixel on map image. I've Manhattan image map rotated 28.34. I could not find correct formula for converting lat/lng to pixel and vice versa. I May be missing something important?

Manhattan map image (1816x8160) having following lat/lon of corners.

TopLeft: (-73.9308,40.8883)
TopRight: (-73.8584,40.858)
BottomLeft: (-74.0665,40.7024)
BottomRight: (-73.9944,40.6718)
Map is not true north and rotated at 28.34, also its UTM Zone 18N (78W to 72W). Here are further details about this map taken from PDF Maps iOS app.

Size (pixels): 1816 x 6160
Pixel Resolution: 3.829 meters
Bounds (pixels): (-1624, -3518) x (7866, 7719)
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32618"]]

How to convert lat/lon to x y and vice versa?

Improve this question
5
  • Here is sample map image link
    – Tariq Mahmood
    Commented Jan 27, 2016 at 10:39
  • Could you apply a rotation (linear transform) on the pixels and then convert? Your scale appears to be small enough that each pixel should represent an equal amount of latitude and an equal amount of longitude as other pixels (although the latitude and longitude amounts themselves will be different)
    – user1462
    Commented Jan 27, 2016 at 15:40
  • @barrycarter thanks for your comments. Yes, I tried with rotation 'rotX = X * cos(angle) - Y * sin(angle); ' 'rotY = X * sin(angle) + Y * cos(angle);' then convert to lat/lng but results are still different. I need a pixel perfect solution.
    – Tariq Mahmood
    Commented Feb 1, 2016 at 6:45
  • Note that your pixel resolution is 3.829 meters, which equates to .0000344610 degrees, but your corner points are given to only 4 post-decimal digits of precision. I think this may be causing the small errors you are seeing.
    – user1462
    Commented Feb 2, 2016 at 12:50
  • Also, above lat/lng of corners were manually taken from Qgis after georeferencing map image. From PDF Maps iOS app Lat/lng of corners for same map are. TopLeft: 40.88834, -73.93069 TopRight: 40.85784,-73.85842 BottomLeft: 40.70223, -74.06649 BottomRight: 40.67182,-73.99436
    – Tariq Mahmood
    Commented Feb 2, 2016 at 13:37

2 Answers 2

Reset to default
1

Just as a quick double check, when you rotated the coordinates, you did "rotate from the center", correct? In other words, your resulting coordinates were:

{x Cos[angle] + (2 cy - y) Sin[angle], y Cos[angle] + x Sin[angle]}

where cy is the center y coordinate (ie, half the y width, or 3080 in your case), and angle is the angle of rotation, correct?

When I did that translation on your sample image and compared it to ImageMagick's rotated version of the image, I got pixel perfect accuracy (after rounding):

enter image description here

(note: I can't see this image myself, but I can when I download it, so I'll assume it's just me; if no one else can see it either, please let me know).

[old "answer" below]

This is probably not helpful, but too long for a comment.

The formulas I found using https://github.com/barrycarter/bcapps/blob/master/STACK/bc-solve-gis-178354.m are:

{-73.9308 + 0.00003988980716253542*x - 0.000016631940188748616*y, 
 40.8883 - 0.000016694214876035256*x - 0.000022784654982228672*y}

to convert pixel x,y to lon,lat, and

{1.9928145063077821*10^6 - 14017.262662576404*lat + 19202.720184031*lon, 
 334427.99241686985 - 33618.80202858345*lat - 14069.747257820272*lon}

to convert lon,lat to x,y.

These are close, but don't quite match what you need.

Having better precision on the corner coordinates may allow refinement of these formulas.

Improve this answer
3
  • Thanks @barrycarter, but as you say results are quite different.
    – Tariq Mahmood
    Commented Feb 2, 2016 at 12:45
  • I've tested it but as you say results are quite different. For example at x=930,y=3590 lat/lng should be lat=40.764296,lng=-73.973027. But using these formulas pixelToLngLat(930,3590)= -73.95341114, 40.79097747 and latLngToPixel(40.764296,-73.973027)= 927.323374, 4762.989144 which are quite different results.
    – Tariq Mahmood
    Commented Feb 2, 2016 at 13:03
  • thanks for pointing to "Rotate from center". I've seen some good results with it. I'll further test and let you know.
    – Tariq Mahmood
    Commented Feb 4, 2016 at 7:10
0

Try these: 

def convert_latlon_to_pixelXY(lat,lon,prj="",jgw=[]):
    dst = osr.SpatialReference(prj)
    transform = jgw
    if isTif:
        x_NorthWest = transform[0]
        y_NorthWest = transform[3]
        meter_perPixel_X = transform[1]
        meter_perPixel_Y = transform[5]
    else:
        x_NorthWest = float(transform[4])
        y_NorthWest = float(transform[5])
        meter_perPixel_X = float(transform[0])
        meter_perPixel_Y = float(transform[3])

src = osr.SpatialReference()
src.SetWellKnownGeogCS('WGS84')
ct = osr.CoordinateTransformation(src,dst)
xy = ct.TransformPoint(lon, lat)
Improve this answer
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.