How to convert lat, long to x y for plotting on the map?

I have an 14080*13624 .jpg image that is the country map.

also I know the lat and long of 4 corner of the map.

I have the lat and long of some other point that I want to scatter them on my map.

I used the code below but the result dosen't match the real position of points.

img = imread('map.jpg');
imshow(img);
hold on
%% Desiger point
lat = 32.778473;
lon = 51.648468;
h = 400;
%% center
lat0 = 31.8978;
lon0 = 53.6133;
h0 = 400;
%% 
wgs84 = wgs84Ellipsoid;
[xEast,yNorth,zUp] = geodetic2enu(lat,lon,h,lat0,lon0,h0,wgs84);
pix_size_x = 14080/2;
pix_size_y = 13824/2;
stepx = (8.2281e+05)/pix_size_x;
stepy = (9.2292e+05)/pix_size_y;
x = pix_size_x + xEast/stepx;
y = pix_size_y - yNorth/stepy;
scatter(x,y)
axis equal

Would you please help me solve this problem?!

[Walter Roberson]

pix_size_x = 14080/2;
pix_size_y = 13824/2;
stepx = (8.2281e+05)/pix_size_x;
stepy = (9.2292e+05)/pix_size_y;

Where are those numbers from? If the pix_size have to do with the image size then use size() to calculate them instead of hardcoding them.

[Sabereh]

@Walter Roberson

8.2281e+5 and 9.2292e+5 are the xEast and yNorth when lat and long equal to one of the corners.

[Walter Roberson]

Cleaning up the code a little but not solving the problem.

It is not immediately obvious to me why you divide the pixels by 2. Something to do with the reference point being at the center of the map?

I would suggest that you have the coordinates of the bounds of the map as input, and let it calculate the centre and span.

%% map size
long_covered = 8.2281e+05;
lat_covered = 9.2292e+05;

%% Desiger point
lat = 32.778473;
lon = 51.648468;
h = 400;
%% center
lat0 = 31.8978;
lon0 = 53.6133;
h0 = 400;

%% 
img = imread('map.jpg');
imshow(img);
hold on

wgs84 = wgs84Ellipsoid;
[xEast,yNorth,zUp] = geodetic2enu(lat,lon,h,lat0,lon0,h0,wgs84);
pix_size_x = size(img,2)/2;
pix_size_y = size(img,1)/2;
stepx = long_covered/pix_size_x;
stepy = lat_covered/pix_size_y;
x = pix_size_x + xEast/stepx;
y = pix_size_y - yNorth/stepy;
scatter(x,y)
axis equal

[Sabereh]

@Walter Roberson yes,the divided by 2 is because of the refrence point is center.

[Benjamin Thompson]

If you have the Mapping Toolbox you can use makerefmat and latlon2pix.