[GIS] convert lat/lng to pixel and vice versa

coordinate systemlatitude longituderesampling

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 (1816×8160) 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?

Best Answer

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.

Related Question