Convert Lat Long to X Y having lat/lon of corners of a map image

Question

I've map image (1816 x 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? Any help will be much appreciated. Thanks in advance.

var dot_size = 15;
var longitude_shift = 0; //-28.34;   // number of pixels your map's prime meridian is off-center.
var x_pos = 0; //54;
var y_pos = 0; //19;
var map_width = 1380; //1816; //430;
var map_height = 4682; //6160; //332;
var half_dot = Math.floor(dot_size / 2);
// Converts from degrees to radians.
Math.radians = function(degrees) {
  return degrees * Math.PI / 180;
}; 
// Converts from radians to degrees.
Math.degrees = function(radians) {
  return radians * 180 / Math.PI;
};

/* N 40.88839 -73.93071 //-73.9308
E 40.85789 -73.85843  //40.858 -73.8584
W 40.70228 -74.06652 //40.7024 -74.0665
S 40.67185 -73.99437 //40.6718 -73.9944 */
var bottomX = 40.67185;
var bottomY = -73.99437;
var topX = 40.88839; //-73.9308; //-73.9308,40.8883
var topY = -73.93071;

var degreesPerPixelX = (bottomX - topX) / map_width; //0.07225 / map_width;
var degreesPerPixelY =  (bottomY - topY) / map_height; //0.18605/ map_height; 

// These should roughly box Germany - use the actual values appropriate to your image
var minLat = bottomX;
var minLong = bottomY;
var maxLat = topX;
var maxLong = topY;

// Map image size (in points)
var mapSize = {'width': map_width, 'height': map_height};

// Determine the map scale (points per degree)
var xScale = mapSize.width / (maxLong - minLong);
var yScale = mapSize.height / (maxLat - minLat);

var south = Math.radians(40.67185); //lat 47.2
var north = Math.radians(40.88839); //lat 55.2
var west = Math.radians(-74.06652); //long 5.8
var east = Math.radians(-73.85843); //long 15.2

// Formula for mercator projection y coordinate:
function mercY(lat) { return Math.log(Math.tan(lat/2 + Math.PI/4)); }

// Some constants to relate chosen area to screen coordinates
var ymin = mercY(south);
var ymax = mercY(north);
var xFactor = mapSize.width/(east - west);
var yFactor = mapSize.height/(ymax - ymin);

var mapLonLeft = -74.06652; //9.8;
var mapLonRight = -73.85843; //10.2;
var mapLonDelta = mapLonRight - mapLonLeft;

mapLatBottom = 40.67185; //53.45;
mapLatBottomRadian = mapLatBottom * Math.PI / 180;

function convertGeoToPixel(lat, lon)
{
    pX = (lon - mapLonLeft) * (mapSize.width / mapLonDelta);
    lat1 = lat * Math.PI / 180;
    worldMapWidth = ((mapSize.width / mapLonDelta) * 360) / (2 * Math.PI);
    mapOffsetY = (worldMapWidth / 2 * Math.log((1 + Math.sin(mapLatBottomRadian)) / (1 - Math.sin(mapLatBottomRadian))));
    pY = mapSize.height - ((worldMapWidth / 2 * Math.log((1 + Math.sin(lat1)) / (1 - Math.sin(lat1)))) - mapOffsetY);

    return 'x:'+pX+', y:'+pY;
}
function convertPixelToGeo(tx, ty)
{   
    /* called worldMapWidth in Raphael's Code, but I think that's the radius since it's the map width or circumference divided by 2*PI  */   
    var worldMapRadius = mapSize.width / mapLonDelta * 360/(2 * Math.PI);     
    var mapOffsetY = ( worldMapRadius / 2 * Math.log( (1 + Math.sin(mapLatBottomRadian) ) / (1 - Math.sin(mapLatBottomRadian))  ));
    var equatorY = mapSize.height + mapOffsetY;   
    var a = (equatorY-ty)/worldMapRadius;

    var lat = 180/Math.PI * (2 * Math.atan(Math.exp(a)) - Math.PI/2);
    var long = mapLonLeft+tx/mapSize.width*mapLonDelta;
    return 'lat:'+lat+', lng:'+long;
}
function draw_point(x, y) {
    dot = '<div style="position:absolute;width:' + dot_size + 'px;height:' + dot_size + 'px;top:' + y + 'px;left:' + x + 'px;background:#00ff00"></div>';
    document.body.innerHTML += dot;
}
function plot_point(lat, lng) {
 spotLat = lat;
 spotLong = lng;
    // Mercator projection
    // longitude: just scale and shift
    x = (map_width * (180 + lng) / 360) % map_width + longitude_shift;

    // latitude: using the Mercator projection
    lat1 = lat * Math.PI / 180;  // convert from degrees to radians
    y = Math.log(Math.tan((lat1/2) + (Math.PI/4)));  // do the Mercator projection (w/ equator of 2pi units)
    y = (map_height / 2) - (map_width * y / (2 * Math.PI)) + y_pos;   // fit it to our map

    x -= x_pos;
    y -= y_pos;
 
 // position of map image for point
 //var newXY = 'x:' (spotLong - minLong) * xScale + ', y:' + (spotLat - minLat) * yScale +'<br/>';
 //var y = (spotLat - minLat) * yScale;
 
 //alert('x: ' + kavraX(Math.radians(lat),Math.radians(lng)) + ', y: ' + kavraY(Math.radians(lat),Math.radians(lng)));
 
 strText = 'kavra x:' + kavraX(Math.radians(lat),Math.radians(lng)) + ', y:' + kavraY(Math.radians(lat),Math.radians(lng)) + '<br/>';
 strText += 'x:' + x + ', y:' + y + '<br/>';
 strText += 'x:'+(spotLong - minLong) * xScale +', y:' + (spotLat - minLat) * yScale +'<br/>';
 strText += 'x:'+((Math.radians(lng) - west)*xFactor)+' ,y:'+((ymax - mercY(Math.radians(lat)))*yFactor)+'<br/>';
 strText += convertGeoToPixel(lat,lng)+'<br/>' ;
 //floatingDiv = '<div style="position:fixed;top:10px;left:10px;">'+strText+'</div>';
 //document.body.innerHTML += floatingDiv;
 $('#leftDiv').html(strText);
    draw_point(x - half_dot, y - half_dot);
}
function kavraX (latitude, longitude) // Kavra for Kavrayskiy 
// formula from http://en.wikipedia.org/wiki/Kavrayskiy_VII_projection
{
  return ((3 * longitude) / 2 * Math.PI)*Math.sqrt(Math.pow(Math.PI, 2)/3 - Math.pow(latitude, 2));
} 
function kavraY (latitude, longitude) 
{
  return latitude*-1;
}

$(document).ready(function() {
//-73.949321, 40.796997
plot_point(40.764296, -73.973027);
  $('img').click(function(e) {
    var offset = $(this).offset();
    var relativeX = (e.pageX - offset.left);
 var relativeY = (e.pageY - offset.top);
 var clickedLon = topX + relativeX * degreesPerPixelX;
    var clickedLat = bottomY + relativeY * degreesPerPixelY;
 alert(relativeX+':'+relativeY+' lat:'+clickedLat+", lon:"+clickedLon);
  });
  
  $('#parentDiv').mousemove(function(e) {
    var offset = $('img').offset();
    var relativeX = (e.pageX - offset.left);
 var relativeY = (e.pageY - offset.top);
 var clickedLat = topX + relativeX * degreesPerPixelX;
    var clickedLon = topY + relativeY * degreesPerPixelY;
 //alert(relativeX+':'+relativeY+' lat:'+clickedLat+", lon:"+clickedLon);
 var strText = relativeX+':'+relativeY+' lat:'+clickedLat+", lon:"+clickedLon+'<br/>';
 strText += 'lat:'+(relativeY / yScale + minLat)+', long:'+(relativeX / xScale + minLong)+'<br/>';
 strText += convertPixelToGeo(relativeX,relativeY)+'<br/>';
 //floatingDiv = '<div style="position:fixed;top:10px;right:10px;">'+strText+'</div>';
 //document.body.innerHTML += floatingDiv;
 $('#rightDiv').html(strText);
  });
});

/*$(function() {
$("#test").click(function(e) {
  var offset = $(this).offset();
  var relativeX = (e.pageX - offset.left);
  var relativeY = (e.pageY - offset.top);
  alert(relativeX+':'+relativeY);
  $(".position").val("afaf");
});
});*/
   

function onClick (evt) {
 alert(evt.pageX +' '+ evt.pageY);
    var x = getEventOffsetFromImageLeft(evt);
    var y = getEventOffsetFromImageTop(evt);

    var clickedLon = topX + x * degreesPerPixelX;
    var clickedLat = bottomY + y * degreesPerPixelY;
}
</script>
</head>

<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script>
<html>
<head>

</head>
<!-- onload="plot_point(40.756, -73.986)" -->
<body >
<div id='parentDiv'>
 <!-- image found at https://i.stack.imgur.com/BXgSw.jpg -->
    <img src="https://i.stack.imgur.com/BXgSw.jpg" style="position:absolute;top:0px;left:0px" >
 <div id="leftDiv" style="position:fixed;top:10px;left:10px;"></div>
 <div id="rightDiv" style="position:fixed;top:10px;right:10px;"></div>
 </div>
</body>
</html>

Take a look at this, http://stackoverflow.com/questions/19849854/get-latitude-and-longitude-from-x-and-y-coordinates — Alvin Varghese, Jan 09 '16 at 16:41
I tried and its giving x: -11.036900291825647, y: -0.012414987560827094 for 40.756, -73.986? — Tariq Mahmood, Jan 09 '16 at 22:32
Run Code Snippet added you can try yourself. Remember source image resolution is changed by i.stack.imggur so map img width/height will change. 'var map_width = 1380; //1816; var map_height = 4682; //6160; ' — Tariq Mahmood, Jan 11 '16 at 11:54
In code snippet above you can see I tried different approaches but none giving me exact results. Anyone can try and compare result with google maps or other maps applications. — Tariq Mahmood, Jan 11 '16 at 12:43
It looks like results are not correct because my map is rotated 28.34. Anyone have thought on how to get x, y from lat,long (and vice versa) on rotated map? — Tariq Mahmood, Jan 11 '16 at 23:23
Re-posted as http://stackoverflow.com/questions/34853937/convert-lat-lng-to-pixel#comment57817373_34853937 — PolyGeo, Jan 27 '16 at 19:47

Convert Lat Long to X Y having lat/lon of corners of a map image

0 Answers0