Get a streetmap of where you live using OpenStreetMap.org (via Google Geocaching)
andyfaff
Sun, 03/08/2009 - 03:52 am
To use type:
getPlotOfAddress("New Illawarra Rd Sydney Australia", 10000)
Where the first string is a search string you would use on googlemaps. The second number is the length of the image square you want.
It works by using Geoff Duttons code to work out the lat/long of the address you enter using Google geocaching. Then this lat-long is sent to Openstreetmap.org to get a maptile. It works fairly well but needs some experimentation with the zoom level and the scale of the map requested.
Please use responsibly, don't overdo it and get blacklisted by Google/OpenStreetmap by trying to programatically download maps of the world.
#pragma rtGlobals=1 // Use modern global access method.
CONSTANT DEG2RAD = 0.01745329;
CONSTANT RAD2DEG = 57.29577951;
Structure GPSpos
variable UTMNorthing
variable UTMEasting
string UTMzone
variable Lat
variable Long
Endstructure
Function getPlotOfAddress(address, zoomlevel)
string address
variable zoomlevel
variable minLong, maxLong, minLat,maxLat
struct GPSpos gps
struct GPSpos mapArea
//get the lat/long of the address
Geocodewithgoogle(address, gps)
//convert to UTM so that we can get a zoom box
LLtoUTM(22, gps.lat, gps.Long, gps)
minLat = gps.UTMNorthing-zoomLevel/2
maxLat = gps.UTMNorthing+zoomLevel/2
minLong = gps.UTMEasting-zoomLevel/2
maxLong = gps.UTMEasting+zoomLevel/2
//convert the min long in UTM to LL
UTMtoLL(22, minLong, minLat, gps.UTMzone, maparea)
minlong = mapArea.long
minlat = mapArea.lat
//convert the max long in UTM to LL
UTMtoLL(22, maxLong, maxLat, gps.UTMzone, maparea)
maxlong = mapArea.long
maxlat = mapArea.lat
string url="", fiLoc=""
fiLoc = SpecialDirpath("Temporary",0,0,0 )+"abcd.png"
sprintf url, "http://www.openstreetmap.org/export/finish?maxlat=%g&minlon=%g&maxlon=%…;,maxLat, minLong, maxLong,minLat
easyHttp/File=filoc url
ImageLoad/T=png/O/G fiLoc
End
Function LLtoUTM( ReferenceEllipsoid, Lat, Long, gps)
//converts lat/long to UTM coords. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees
//Written by Chuck Gantz- chuck.gantz@globalstar.com
variable ReferenceEllipsoid, Lat, Long
struct GPSpos &gps
variable UTMEasting, UTMNorthing
string UTMZone
make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
make/o/n=23/T EllipsoidName
EquatorialRadius={6377563, 6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough", "international", "Krassovsky", "Modified Airy", "Modified Everest", "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}
variable a = EquatorialRadius[ReferenceEllipsoid]
variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
variable k0 = 0.9996;
variable LongOrigin;
variable eccPrimeSquared;
variable N, T, C, AA, M;
//Make sure the longitude is between -180.00 .. 179.9
variable LongTemp = (Long+180)-floor((Long+180)/360)*360-180; // -180.00 .. 179.9;
variable LatRad = Lat*DEG2RAD;
variable LongRad = LongTemp*DEG2RAD;
variable LongOriginRad;
variable ZoneNumber;
ZoneNumber = floor((LongTemp + 180)/6) + 1;
if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
ZoneNumber = 32;
endif
// Special zones for Svalbard
if( Lat >= 72.0 && Lat < 84.0 )
if ( LongTemp >= 0.0 && LongTemp < 9.0 )
ZoneNumber = 31;
elseif( LongTemp >= 9.0 && LongTemp < 21.0 )
ZoneNumber = 33;
elseif(LongTemp >= 21.0 && LongTemp < 33.0 )
ZoneNumber = 35;
elseif(LongTemp >= 33.0 && LongTemp < 42.0 )
ZoneNumber = 37;
endif
endif
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
LongOriginRad = LongOrigin * DEG2RAD;
//compute the UTM Zone from the latitude and longitude
eccPrimeSquared = (eccSquared)/(1-eccSquared);
N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
T = tan(LatRad)*tan(LatRad);
C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
AA = cos(LatRad)*(LongRad-LongOriginRad);
M = (1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
M -= (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
M += (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
M -= (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)
M *= a
UTMEasting = (k0*N*(AA+(1-T+C)*AA*AA*AA/6+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*AA*AA*AA*AA*AA/120)+ 500000.0);
UTMNorthing = (k0*(M+N*tan(LatRad)*(AA*AA/2+(5-T+9*C+4*C*C)*AA*AA*AA*AA/24+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*AA*AA*AA*AA*AA*AA/720)));
if(Lat < 0)
UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
endif
gps.UTMEasting = UTMEasting
gps.UTMNorthing = UTMNorthing
gps.UTMZone = num2str(ZoneNumber)+ UTMLetterDesignator(Lat)
End
Function/t UTMLetterDesignator( Lat)
variable Lat
//This routine determines the correct UTM letter designator for the given latitude
//returns "Z" if latitude is outside the UTM limits of 84N to 80S
//Written by Chuck Gantz- chuck.gantz@globalstar.com
string LetterDesignator;
if((84 >= Lat) && (Lat >= 72))
LetterDesignator = "X";
elseif((72 > Lat) && (Lat >= 64))
LetterDesignator = "W";
elseif((64 > Lat) && (Lat >= 56))
LetterDesignator = "V";
elseif((56 > Lat) && (Lat >= 48))
LetterDesignator = "U";
elseif((48 > Lat) && (Lat >= 40))
LetterDesignator = "T";
elseif((40 > Lat) && (Lat >= 32))
LetterDesignator = "S";
elseif((32 > Lat) && (Lat >= 24))
LetterDesignator = "R";
elseif((24 > Lat) && (Lat >= 16))
LetterDesignator = "Q";
elseif((16 > Lat) && (Lat >= 8))
LetterDesignator = "P";
elseif(( 8 > Lat) && (Lat >= 0))
LetterDesignator = "N";
elseif(( 0 > Lat) && (Lat >= -8))
LetterDesignator = "M";
elseif((-8> Lat) && (Lat >= -16))
LetterDesignator = "L";
elseif((-16 > Lat) && (Lat >= -24))
LetterDesignator = "K";
elseif((-24 > Lat) && (Lat >= -32))
LetterDesignator = "J";
elseif((-32 > Lat) && (Lat >= -40))
LetterDesignator = "H";
elseif((-40 > Lat) && (Lat >= -48))
LetterDesignator = "G";
elseif((-48 > Lat) && (Lat >= -56))
LetterDesignator = "F";
elseif((-56 > Lat) && (Lat >= -64))
LetterDesignator = "E";
elseif((-64 > Lat) && (Lat >= -72))
LetterDesignator = "D";
elseif((-72 > Lat) && (Lat >= -80))
LetterDesignator = "C";
else
LetterDesignator = "Z"; //This is here as an error flag to show that the Latitude is outside the UTM limits
endif
return LetterDesignator;
End
Function UTMtoLL(ReferenceEllipsoid, UTMEasting, UTMNorthing, UTMZone, gps)
variable referenceEllipsoid,UTMeasting,UTMNorthing
string UTMzone
struct GPSpos &gps
Variable Lat, Long
//converts UTM coords to lat/long. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees.
//Written by Chuck Gantz- chuck.gantz@globalstar.com
make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
make/o/n=23/T EllipsoidName
EquatorialRadius={6377563, 6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough", "international", "Krassovsky", "Modified Airy", "Modified Everest", "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}
variable k0 = 0.9996;
variable a = EquatorialRadius[ReferenceEllipsoid];
variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
variable eccPrimeSquared;
variable e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
variable N1, T1, C1, R1, D, M;
variable LongOrigin;
variable mu, phi1, phi1Rad;
variable x, y;
variable ZoneNumber;
string ZoneLetter;
variable NorthernHemisphere; //1 for northern hemispher, 0 for southern
x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
y = UTMNorthing;
//UTMZone is something like 32P. Zone number is 32, Zone letter is P
sscanf UTMZone, "%d%s", ZoneNumber, ZoneLetter
if((char2num(ZoneLetter) - char2num("N")) >= 0)
NorthernHemisphere = 1;//povariable is in northern hemisphere
else
NorthernHemisphere = 0;//povariable is in southern hemisphere
y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
endif
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
eccPrimeSquared = (eccSquared)/(1-eccSquared);
M = y / k0;
mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)+ (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)+(151*e1*e1*e1/96)*sin(6*mu);
phi1 = phi1Rad*rad2deg;
N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
T1 = tan(phi1Rad)*tan(phi1Rad);
C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
R1 = a*(1-eccSquared)/((1-eccSquared*sin(phi1Rad)*sin(phi1Rad))^1.5);
D = x/(N1*k0);
Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
Lat = Lat * rad2deg;
Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)*D*D*D*D*D/120)/cos(phi1Rad);
Long = LongOrigin + Long * rad2deg;
gps.Long = Long
gps.Lat = Lat
End
//Function Uses Google maps to return latitude and longitude of an address
// Requires: easyHttp XOP
// Geoff Dutton March 5, 2009
Function GeocodeWithGoogle(address, gps)
string address
struct GPSpos &gps
variable lat, long
String url, GmapOutput
// build a URL for Google maps
address = urlencode(Address)
Sprintf url, "http://maps.google.com/maps/geo?q=%s&output=csv", address
// use easyHttp to contact Google maps, Google maps returns lat, long encoded text
easyHttp url
GmapOutput = S_getHttp
// create a list from returned Google text
GmapOutput = ReplaceString(",", GmapOutput, ";")
// Google returns 200 if address was found
if ( str2num(StringFromList( 0, GmapOutput)) == 200 )
lat = str2num(StringFromList(2, GmapOutput))
long = str2num(StringFromList(3, GmapOutput))
endif
gps.Lat = lat
gps.long = long
end
CONSTANT DEG2RAD = 0.01745329;
CONSTANT RAD2DEG = 57.29577951;
Structure GPSpos
variable UTMNorthing
variable UTMEasting
string UTMzone
variable Lat
variable Long
Endstructure
Function getPlotOfAddress(address, zoomlevel)
string address
variable zoomlevel
variable minLong, maxLong, minLat,maxLat
struct GPSpos gps
struct GPSpos mapArea
//get the lat/long of the address
Geocodewithgoogle(address, gps)
//convert to UTM so that we can get a zoom box
LLtoUTM(22, gps.lat, gps.Long, gps)
minLat = gps.UTMNorthing-zoomLevel/2
maxLat = gps.UTMNorthing+zoomLevel/2
minLong = gps.UTMEasting-zoomLevel/2
maxLong = gps.UTMEasting+zoomLevel/2
//convert the min long in UTM to LL
UTMtoLL(22, minLong, minLat, gps.UTMzone, maparea)
minlong = mapArea.long
minlat = mapArea.lat
//convert the max long in UTM to LL
UTMtoLL(22, maxLong, maxLat, gps.UTMzone, maparea)
maxlong = mapArea.long
maxlat = mapArea.lat
string url="", fiLoc=""
fiLoc = SpecialDirpath("Temporary",0,0,0 )+"abcd.png"
sprintf url, "http://www.openstreetmap.org/export/finish?maxlat=%g&minlon=%g&maxlon=%…;,maxLat, minLong, maxLong,minLat
easyHttp/File=filoc url
ImageLoad/T=png/O/G fiLoc
End
Function LLtoUTM( ReferenceEllipsoid, Lat, Long, gps)
//converts lat/long to UTM coords. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees
//Written by Chuck Gantz- chuck.gantz@globalstar.com
variable ReferenceEllipsoid, Lat, Long
struct GPSpos &gps
variable UTMEasting, UTMNorthing
string UTMZone
make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
make/o/n=23/T EllipsoidName
EquatorialRadius={6377563, 6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough", "international", "Krassovsky", "Modified Airy", "Modified Everest", "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}
variable a = EquatorialRadius[ReferenceEllipsoid]
variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
variable k0 = 0.9996;
variable LongOrigin;
variable eccPrimeSquared;
variable N, T, C, AA, M;
//Make sure the longitude is between -180.00 .. 179.9
variable LongTemp = (Long+180)-floor((Long+180)/360)*360-180; // -180.00 .. 179.9;
variable LatRad = Lat*DEG2RAD;
variable LongRad = LongTemp*DEG2RAD;
variable LongOriginRad;
variable ZoneNumber;
ZoneNumber = floor((LongTemp + 180)/6) + 1;
if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
ZoneNumber = 32;
endif
// Special zones for Svalbard
if( Lat >= 72.0 && Lat < 84.0 )
if ( LongTemp >= 0.0 && LongTemp < 9.0 )
ZoneNumber = 31;
elseif( LongTemp >= 9.0 && LongTemp < 21.0 )
ZoneNumber = 33;
elseif(LongTemp >= 21.0 && LongTemp < 33.0 )
ZoneNumber = 35;
elseif(LongTemp >= 33.0 && LongTemp < 42.0 )
ZoneNumber = 37;
endif
endif
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
LongOriginRad = LongOrigin * DEG2RAD;
//compute the UTM Zone from the latitude and longitude
eccPrimeSquared = (eccSquared)/(1-eccSquared);
N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
T = tan(LatRad)*tan(LatRad);
C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
AA = cos(LatRad)*(LongRad-LongOriginRad);
M = (1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
M -= (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
M += (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
M -= (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)
M *= a
UTMEasting = (k0*N*(AA+(1-T+C)*AA*AA*AA/6+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*AA*AA*AA*AA*AA/120)+ 500000.0);
UTMNorthing = (k0*(M+N*tan(LatRad)*(AA*AA/2+(5-T+9*C+4*C*C)*AA*AA*AA*AA/24+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*AA*AA*AA*AA*AA*AA/720)));
if(Lat < 0)
UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
endif
gps.UTMEasting = UTMEasting
gps.UTMNorthing = UTMNorthing
gps.UTMZone = num2str(ZoneNumber)+ UTMLetterDesignator(Lat)
End
Function/t UTMLetterDesignator( Lat)
variable Lat
//This routine determines the correct UTM letter designator for the given latitude
//returns "Z" if latitude is outside the UTM limits of 84N to 80S
//Written by Chuck Gantz- chuck.gantz@globalstar.com
string LetterDesignator;
if((84 >= Lat) && (Lat >= 72))
LetterDesignator = "X";
elseif((72 > Lat) && (Lat >= 64))
LetterDesignator = "W";
elseif((64 > Lat) && (Lat >= 56))
LetterDesignator = "V";
elseif((56 > Lat) && (Lat >= 48))
LetterDesignator = "U";
elseif((48 > Lat) && (Lat >= 40))
LetterDesignator = "T";
elseif((40 > Lat) && (Lat >= 32))
LetterDesignator = "S";
elseif((32 > Lat) && (Lat >= 24))
LetterDesignator = "R";
elseif((24 > Lat) && (Lat >= 16))
LetterDesignator = "Q";
elseif((16 > Lat) && (Lat >= 8))
LetterDesignator = "P";
elseif(( 8 > Lat) && (Lat >= 0))
LetterDesignator = "N";
elseif(( 0 > Lat) && (Lat >= -8))
LetterDesignator = "M";
elseif((-8> Lat) && (Lat >= -16))
LetterDesignator = "L";
elseif((-16 > Lat) && (Lat >= -24))
LetterDesignator = "K";
elseif((-24 > Lat) && (Lat >= -32))
LetterDesignator = "J";
elseif((-32 > Lat) && (Lat >= -40))
LetterDesignator = "H";
elseif((-40 > Lat) && (Lat >= -48))
LetterDesignator = "G";
elseif((-48 > Lat) && (Lat >= -56))
LetterDesignator = "F";
elseif((-56 > Lat) && (Lat >= -64))
LetterDesignator = "E";
elseif((-64 > Lat) && (Lat >= -72))
LetterDesignator = "D";
elseif((-72 > Lat) && (Lat >= -80))
LetterDesignator = "C";
else
LetterDesignator = "Z"; //This is here as an error flag to show that the Latitude is outside the UTM limits
endif
return LetterDesignator;
End
Function UTMtoLL(ReferenceEllipsoid, UTMEasting, UTMNorthing, UTMZone, gps)
variable referenceEllipsoid,UTMeasting,UTMNorthing
string UTMzone
struct GPSpos &gps
Variable Lat, Long
//converts UTM coords to lat/long. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees.
//Written by Chuck Gantz- chuck.gantz@globalstar.com
make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
make/o/n=23/T EllipsoidName
EquatorialRadius={6377563, 6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough", "international", "Krassovsky", "Modified Airy", "Modified Everest", "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}
variable k0 = 0.9996;
variable a = EquatorialRadius[ReferenceEllipsoid];
variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
variable eccPrimeSquared;
variable e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
variable N1, T1, C1, R1, D, M;
variable LongOrigin;
variable mu, phi1, phi1Rad;
variable x, y;
variable ZoneNumber;
string ZoneLetter;
variable NorthernHemisphere; //1 for northern hemispher, 0 for southern
x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
y = UTMNorthing;
//UTMZone is something like 32P. Zone number is 32, Zone letter is P
sscanf UTMZone, "%d%s", ZoneNumber, ZoneLetter
if((char2num(ZoneLetter) - char2num("N")) >= 0)
NorthernHemisphere = 1;//povariable is in northern hemisphere
else
NorthernHemisphere = 0;//povariable is in southern hemisphere
y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
endif
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
eccPrimeSquared = (eccSquared)/(1-eccSquared);
M = y / k0;
mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)+ (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)+(151*e1*e1*e1/96)*sin(6*mu);
phi1 = phi1Rad*rad2deg;
N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
T1 = tan(phi1Rad)*tan(phi1Rad);
C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
R1 = a*(1-eccSquared)/((1-eccSquared*sin(phi1Rad)*sin(phi1Rad))^1.5);
D = x/(N1*k0);
Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
Lat = Lat * rad2deg;
Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)*D*D*D*D*D/120)/cos(phi1Rad);
Long = LongOrigin + Long * rad2deg;
gps.Long = Long
gps.Lat = Lat
End
//Function Uses Google maps to return latitude and longitude of an address
// Requires: easyHttp XOP
// Geoff Dutton March 5, 2009
Function GeocodeWithGoogle(address, gps)
string address
struct GPSpos &gps
variable lat, long
String url, GmapOutput
// build a URL for Google maps
address = urlencode(Address)
Sprintf url, "http://maps.google.com/maps/geo?q=%s&output=csv", address
// use easyHttp to contact Google maps, Google maps returns lat, long encoded text
easyHttp url
GmapOutput = S_getHttp
// create a list from returned Google text
GmapOutput = ReplaceString(",", GmapOutput, ";")
// Google returns 200 if address was found
if ( str2num(StringFromList( 0, GmapOutput)) == 200 )
lat = str2num(StringFromList(2, GmapOutput))
long = str2num(StringFromList(3, GmapOutput))
endif
gps.Lat = lat
gps.long = long
end
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More