Mega Code Archive
Albers Equal Area formula
<![CDATA[
Title: Albers Equal-Area formula
Question: Formula for converting lat/long data to Albers Equal-Area values to draw maps of an area. The reverse formula to display the location of the mouse being moved around on the map is also provided.
Answer:
{
This was written with Delphi 6, but should work with little or no changes for Delphi 4 or 5 also. Open a new application, put a Button, ScrollBox, Image (in the ScrollBox), and a Label on the form. Cut and paste this listed code and it should compile and run. If there is interest, I'll send in other projections also.
You don't have to know that much about the math, just send the raw latitude/longitude, in degrees and fraction of degrees, to the procedure and mark or draw the line to the returned screen location. Fraction of degrees would be the degrees + 60 / minutes + 3600 / seconds. Should be pretty simple for most people to draw the maps they desire. You will have to get the lat/long files to use from the internet. I'd have sent them also, but they are quite large for this forum.
}
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls, Math;
type
TForm1 = class(TForm)
Button1 : TButton;
ScrollBox1 : TScrollBox;
Image1 : TImage;
Label1 : TLabel;
procedure Button1Click(Sender : TObject);
procedure FormCreate(Sender : TObject);
procedure Image1MouseMove(Sender : TObject; Shift : TShiftState;
X,Y : Integer);
private
{ Private declarations }
public
{ Public declarations }
end;
const
MidLat = 23.0; // middle latitude of the projection.
MidLon = -96.0; // middle longitude of the projection.
SP1 = 29.5; // standard parallel 1.
SP2 = 45.5; // standard parallel 2.
R = 1000.0; // corresponding scale of map.
var
Form1 : TForm1;
OffSet : Integer;
implementation
// These two formulas should give the ability to read in a
// latitude/longitude file and
// display the data in Albers Conical Equal-Area format.
// You can change the scale by changeing the value of 'R'. The
// screen map can be user
// set also. The standard parallels can also be user changed along // with the middle
// latitude/longitude values for your map. And 'OF COURSE' this is
// offered without
// warrenty of anykind. The user assumes all risks using this code
// in their programs.
{$R *.dfm}
// =====================================================================================
// Set screen size.
procedure TForm1.FormCreate(Sender: TObject);
begin
Form1.Image1.Width := 1000;
Form1.Image1.Height := 500;
OffSet := Round(Form1.Image1.Width / 2);
Form1.Caption := 'Albers Conical Equal-Area formulas.';
end;
// =====================================================================================
// X = longitude of point to convert.
// Y = latitude of point to convert.
// MapX = converted longitude location to display on map.
// MapY = converted latitude location to display on map.
// Converts the latitude/longitude to the screen location in Albers Conical Equal-Area.
procedure AlbersEqualAreaSphere(X,Y : Single; var MapX,MapY : Single);
var p : Single; // radius of Latitude circle on conic projection.
Temp1 : Single;
Temp2 : Single;
Temp3 : Single;
n : Single; // ratio of the angle between meridians to the true angle.
begin
n := (Sin(SP1 / 180 * pi) + Sin(SP2 / 180 * pi)) / 2;
Temp1 := Sqr(Cos(SP1 / 180 * pi)) + 2 * n * Sin(SP1 / 180 * pi);
Temp2 := R * Power(Temp1 - 2 * n * Sin(MidLat / 180 * pi),0.5) / n;
Temp3 := n * (X - MidLon);
p := R * Power(Temp1 - 2 * n * Sin(Y / 180 * pi),0.5) / n;
MapX := p * Sin(Temp3 / 180 * pi);
MapY := Temp2 - p * Cos(Temp3 / 180 * pi);
MapX := OffSet + MapX;
MapY := OffSet - MapY;
end;
// =====================================================================================
// Display a arc of circles at 49 degrees north and at 25 degrees
// north, which the
// United States falls within. Also display a few major cities for
// test purposes.
procedure TForm1.Button1Click(Sender: TObject);
var Lp : Integer;
Lon : Single;
X : Single;
Y : Single;
begin
Form1.Image1.Canvas.Pen.Color := clRed;
Lon := -60.0;
for Lp := 1 to 13 do //display arcs, 49 degrees and 25 degrees north,
//65-125 degrees west.
begin
Lon := Lon - 5.0;
AlbersEqualAreaSphere(Lon,49.0,X,Y);
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
AlbersEqualAreaSphere(Lon,25.0,X,Y);
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
end;
// display the locations of 5 cities, for testing of mousemove
//(reverse) formula.
Form1.Image1.Canvas.Pen.Color := clBlack;
AlbersEqualAreaSphere(-73.98,40.77,X,Y); // New York City
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
AlbersEqualAreaSphere(-122.68,37.75,X,Y); // San Francisco
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
AlbersEqualAreaSphere(-104.87,39.75,X,Y); // Denver
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
AlbersEqualAreaSphere(-80.28,25.82,X,Y); // Miami Intl Airport
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
AlbersEqualAreaSphere(-95.35,29.97,X,Y); // Houston
Form1.Image1.Canvas.Ellipse(Round(X-3),Round(Y-3),Round(X+3),Round(Y+3));
end;
// =====================================================================================
// Display the latitude/longitude of the mouse location on the
// screen.
// A reverse for the Albers formula to put the objects on the screen.
procedure TForm1.Image1MouseMove(Sender : TObject; Shift : TShiftState; X,Y : Integer);
var Lat : Single;
Lon : Single;
Test : Single;
Temp1 : Single;
Temp2 : Single;
Temp3 : Single;
n : Single;
p : Single;
SX : Single;
SY : Single;
begin
SX := X - OffSet;
SY := OffSet - Y;
n := (Sin(SP1 / 180 * pi) + Sin(SP2 / 180 * pi)) / 2;
Temp1 := (Cos(SP1 / 180 * pi) * Cos(SP1 / 180 * pi)) + 2 * n * Sin(SP1 / 180 * pi);
Temp2 := R * Power(Temp1 - 2 * n * Sin(MidLat / 180 * pi),0.5) / n;
p := Power((SX * SX) + (Temp2 - SY) * (Temp2 - SY),0.5);
Temp3 := ArcTan(SX / (Temp2 - SY)) * 180 / pi;
Test := (Temp1 - ((p * n / R) * (p * n / R))) / (2 * n);
if (Test 1.0) or (Test Lat := 0
else
Lat := ArcSin((Temp1 - ((p*n/R) * (p*n/R))) / (2*n)) * 180 / pi;
Lon := Temp3 / n + MidLon;
Form1.Label1.Caption := FloatToStr(Lat) + ' ' + FloatToStr(Lon);
end;
// =====================================================================================
end.
]]>