Get a streetmap of where you live using OpenStreetMap.org (via Google Geocaching)

#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=%g&minlat=%g&format=mapnik&mapnik_format=png&mapnik_scale=25000&osmarender_format=png&osmarender_zoom=4&commit=Export",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