Page 1 of 1

Mercator projection

Posted: Sat May 16, 2015 3:19 am
by metalos
Hello, for one of my projects I need a map projected on a point that has Decimal geometric coordinates. The card that I use has a Mercator projection.

Here is the card used:

Image

Here the limits of the map:

North: 52
West: -5.5
East: 10
South: 40.9166

There a few years I found on a forum such procedures:

Code: Select all

Procedure.d mercator_proj(b.d)
  ProcedureReturn(Log(Tan(#PI/4.0+b*#PI/360.0)))
EndProcedure

Macro rad(_deg_)
  ((_deg_)*0.017453292519943295) ; _deg_ * #PI / 180
EndMacro

Procedure.d scale_factor()
  ProcedureReturn(IMAGE_XSIZE/(rad(image_east_longitude)-rad(image_west_longitude)))
EndProcedure

Procedure.d mercator_y(b.d) ; using Latitude
  ProcedureReturn((scale_factor()*(mercator_proj(image_north_latitude)-mercator_proj(b))))
EndProcedure

Procedure.d mercator_x(l.d) ; using Longitude
  ProcedureReturn((scale_factor()*(rad(l)-rad(image_west_longitude))))
EndProcedure
Here are the coordinates of the point I cerche to project on the map:
Latitude: 48.6833
Longitude: 6.2

Only I do not know how to use its procedures to project a point on the map, which is why I need your help. Thank you in advance.

Re: Mercator projection

Posted: Sat May 16, 2015 2:15 pm
by infratec
Since no one answers you:

Code: Select all

EnableExplicit

Procedure.d MercatorProjection(Latitude.d)
  ProcedureReturn Log(Tan(#PI / 4.0 + Latitude * #PI / 360.0))
EndProcedure

Procedure.d MercatorScaleFactor(ImgWidth.i, ImgLongitudeEast.d, ImgLongitudeWest.d)
  ProcedureReturn ImgWidth / (Radian(ImgLongitudeEast) - Radian(ImgLongitudeWest))
EndProcedure

Procedure.d MercatorX(Longitude.d, ImgLongitudeWest.d, ScaleFactor.d)
  ProcedureReturn (Radian(Longitude) - Radian(ImgLongitudeWest)) * ScaleFactor
EndProcedure

Procedure.d MercatorY(Latitude.d, ImgLatitudeNorth.d, ScaleFactor.d)
  ProcedureReturn (MercatorProjection(ImgLatitudeNorth) - MercatorProjection(Latitude)) * ScaleFactor
EndProcedure


Define.i Img, ImgWidth, ImgHeight, X, Y
Define.d Latitude, Longitude, Scale, ImgNorthLatitude, ImgEastLongitude, ImgSouthLatitude, ImgWestLongitude

Latitude = 48.6833
Longitude = 6.2

;Latitude = 49.723359
;Longitude = -1.941223

UsePNGImageDecoder()

Img = CatchImage(#PB_Any, ?MercatorMapImageStart, ?MercatorMapImageEnd - ?MercatorMapImageStart)
If Img
  
  Restore MercatorMapImage
  Read.d ImgNorthLatitude
  Read.d ImgEastLongitude
  Read.d ImgSouthLatitude
  Read.d ImgWestLongitude
  
  ImgWidth = ImageWidth(Img)
  ImgHeight = ImageHeight(Img)
  
  Scale = MercatorScaleFactor(ImgWidth, ImgEastLongitude, ImgWestLongitude)
  X = MercatorX(Longitude, ImgWestLongitude, Scale)
  Y = MercatorY(Latitude, ImgNorthLatitude, Scale)
  
  StartDrawing(ImageOutput(Img))
  Circle(X, Y, 5, $0000FF)
  StopDrawing()
  
  OpenWindow(0, 0, 0, ImgWidth, ImgHeight, "Mercator", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
  ImageGadget(0, 0, 0, 0, 0, ImageID(Img))
  
  Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow
  
EndIf

DataSection
  MercatorMapImage:
  Data.d 52.0, 10.0, 40.9166, -5.5
  MercatorMapImageStart:
  IncludeBinary "MercatorMapImage.png"
  MercatorMapImageEnd:
EndDataSection
I modifed the procedures a bit.
Now they use parameters instead of constants.

Bernd

Re: Mercator projection

Posted: Sat May 16, 2015 3:36 pm
by metalos
It was so simple ... Thank you I can finally forward on my project. Thank you.