[Windows Mobile .NET CF] 找找看附近的停車場資訊 - Day7

[Windows Mobile .NET CF] 找找看附近的停車場資訊 - Day7

看樣子 Hello 開頭的分數不高, 所以要提高檔次了.

取得 GPS 資訊以後, 就可以開始玩 Location based service , 這部分其實台北市政府一直都有在做,
可惜的是, 我相當地懷疑這些政府官員是否有用手機實際用看看這些服務呢?

比方說 , 停車管理工程處就有網頁 http://163.29.129.99/pweb/ ,
台北市政府的即時交通資訊網 : http://its.taipei.gov.tw/atis_index/mapviewer.aspx?Lang=CHT 
上面的網站, 用 PC 上的 IE 看都很辛苦了, 就別提手機了.
當然市政府也有針對手機設計的網頁 : http://its.taipei.gov.tw/m/
不過, 網頁當然沒有 GPS , 沒辦法"找附近的停車場資訊呀”

嗯, 有網頁就可以自己寫程式去讀取資料囉, 所以 HttpWebRequest 實在是個好東西.

在讀取資料之前, 有一個比較困難的地方, 就是台灣的地圖資訊可不是簡單應付的,
以下很多資料參考自 http://wiki.osgeo.org/index.php?title=Taiwan_datums&uselang=zh-tw ,
很大一部分是積丹尼 先生寫的, 想要寫 GPS 程式的人, 強烈推薦這個老外寫的東西. 
真的很令人感動!

另外主要轉換程式也參考自 OGL Library http://ogl-lib.sourceforge.net/index.html

看完這些東西, 應該就知道, 我們從 GPSID 取到的經緯度稱之為 WGS84 的經緯度,
台灣地圖很多是採用 TWD67 或是 TWD97 的經緯度座標以及二度分帶座標, 需要轉換.
之所以會有台灣大地基準以及二度分帶, 三度分帶等這類資訊出現,
不外乎地球是圓的, 地圖卻是平的,
所以我們需要透過麥卡托投影技術, 橢圓體函式 (為什麼國小教科書說地球是圓的?!),
三角函式等等來轉換座標系統.

為了將 WGS84 的經緯度轉換到 TWD97 的二度分帶座標系統, 主要做下面五個步驟:
(這一切似乎是因為我只有找到給 TWD67 用的麥卡扥圓柱投影公式啊?!
公式不難找, 更難找的是公式所用的參數…
不過這樣也好, 不管要 TWD67 還是 TWD97 甚麼座標, 下面都有轉換方法了…)

1. WGS84 轉換成笛卡兒座標系統 , 因為地球是橢圓體, 所以要參考 地球半徑 的算法.
求取卯酉圈之半徑的公式是
N=N(\phi)=\frac{a^2}{\sqrt{(a\cos(\phi))^2+(b\sin(\phi))^2}};\,\! 

以下是程式:
 


private const double WGS84a = 6378137.0d; //the lengths of the semi-major axes of the spheroid
private const double WGS84b = 6356752.314d; //the lengths of the semi-minor axes of the spheroid
private static double WGS84a2 = Math.Pow(WGS84a, 2);
private static double WGS84b2 = Math.Pow(WGS84b, 2);

private const double DegreeToRadianRate = Math.PI / 180.0;

private static double ReturnRadian(double angleValue)
{
    return angleValue * DegreeToRadianRate;
}

/// <summary>
/// WGS84 橢圓體座標轉換笛卡兒座標系統
/// </summary>
/// <param name="Long">Long WGS84 longitude</param>
/// <param name="Lat">Lat WGS84 latitude</param>
/// <param name="h">h WGS84 Antenna altitude above/below mean sea level</param>
/// <returns></returns>
public static DatumPos WGS84toCartesian(double Long, double Lat, double h)
{
    //卯酉圈之半徑 (radius of curvature in prime vertical)
    //N is the radius of curvature in the prime vertical at the point
    //N = a^2 / (a^2 * cos^2(latitude) + b^2 * sin^2(latitude))*(1/2)           
    double theta = ReturnRadian(Lat);

    double N = WGS84a2 /
        Math.Sqrt(Math.Pow(WGS84a * Math.Cos(theta), 2) + Math.Pow(WGS84b * Math.Sin(theta), 2));
    
    //X = (N + h) * cos(1) * cos(2)
    //Y = (N + h) * cos(1) * sin(2)
    //Z = ((b^2/a^2) * N + h) * sin(1)
    //1, 2 are the latitude, longitude of the point 
    //h Height of geoid (mean sea level) above WGS84 Ellipsoidal
    double RadianLong = ReturnRadian(Long);
    double dx = (N + h) * Math.Cos(theta) * Math.Cos(RadianLong);
    double dy = (N + h) * Math.Cos(theta) * Math.Sin(RadianLong);
    double dz = ((WGS84b2 / WGS84a2) * N + h) * Math.Sin(theta);
    return new DatumPos()
    {
        X = dx,
        Y = dy,
        Z = dz
    };
}

2. Bursa-Wolf七參數轉換


/// 七參數 BursaWolf 模型轉換 WGS84 笛卡兒座標 -> TWD67 笛卡兒座標
/// </summary>
/// <param name="wgs84pos"></param>
/// <returns></returns>
public static DatumPos SevenParameter_BursaWolf_WGS84_To_TWD67(DatumPos wgs84pos)
{
    //X67 = △X + SC * (X84 + Y84*Rz - Z84*Ry)
    //Y67 = △Y + SC * (-X84*Rz + Y84 + Z84*Rx)
    //Z67 = △Z + SC * (X84*Ry - Y84*Rx + Z84)
    //RX, RY, RZ = The rotations around the three coordinate axes 
    //SC = The scale difference between the coordinate systems

    /*X67 = DeltaX + SC * (X84 + Y84*Rz - Z84*Ry);
    Y67 = DeltaY + SC * (-X84*Rz + Y84 + Z84*Rx);
    Z67 = DeltaZ + SC * (X84*Ry - Y84*Rx + Z84);*/

    return new DatumPos()
    {
        X = wgs84pos.X + 764.558,
        Y = wgs84pos.Y + 361.299,
        Z = wgs84pos.Z + 178.374
    };
}

3. TWD67 笛卡兒座標轉換橢球體座標系統


/// <summary>
/// Radian convert Angle
/// </summary>
/// <param name="radianValue"></param>
/// <returns></returns>
private static double ReturnAngle(double radianValue)
{
    return radianValue * RadianRateToDegree;
}


private const double TWD67a = 6378160.0d;
private const double TWD67b = 6356774.719d;
private const double m_TWD67aSquared = (TWD67a * TWD67a);
private const double m_TWD67bSquared = (TWD67b * TWD67b);
private const double m_WGS84aSquared = (WGS84a * WGS84a);
private const double m_WGS84bSquared = (WGS84b * WGS84b);
//e^2 = (a^2 - b^2) / a^2
//(e')^2 = e^2 / (1 - e^2)
private const double eccTWD67Squared = (m_TWD67aSquared - m_TWD67bSquared) / m_TWD67aSquared;
private const double eccPrimeTWD67Squared = (eccTWD67Squared) / (1 - eccTWD67Squared);
private const double eccWGS84Squared = ((m_WGS84aSquared - m_WGS84bSquared) / m_WGS84aSquared);
private const double eccPrimeWGS84Squared = (eccWGS84Squared) / (1 - eccWGS84Squared);

/// <summary>
/// TWD67 笛卡兒座標轉換橢球體座標系統
/// </summary>
/// <param name="TWD67"></param>
/// <returns></returns>
public static DatumPos TWD67toEllipsoidal(DatumPos TWD67)
{

    double P = Math.Sqrt(TWD67.X * TWD67.X + TWD67.Y * TWD67.Y);

    double Long67E = ReturnAngle(Math.Atan2(TWD67.Y, TWD67.X));
    double O = Math.Atan2(TWD67.Z * TWD67a, P * TWD67b);
    double Lat67E = ReturnAngle(Math.Atan2(TWD67.Z + eccPrimeTWD67Squared * TWD67b * Math.Pow(Math.Sin(O), 3), P - eccTWD67Squared * TWD67a * Math.Pow(Math.Cos(O), 3)));
    double N = Math.Sqrt(TWD67a / (1 - eccTWD67Squared * Math.Pow(Math.Sin(ReturnRadian(Lat67E)), 2)));
    double h67E = (P / Math.Cos(ReturnRadian(Lat67E))) - N;

    return new DatumPos()
    {
        X = Long67E,
        Y = Lat67E,
        Z = h67E,
    };
}

 

4. 利用麥卡托投影公式, 轉換成二度分帶座標


private const double OffsetX = 250000.0d;    //東距(X軸)    

/// <summary>
/// TWD67(Taiwan Datum 67, TWD-67) Ellipsoidal -> TM2 (N,E)
///UTM: Univeral Transverse Mercator, 橫麥卡托投影經差二度分帶.
///麥卡扥圓柱投影為心射投影, 所得經緯線為正方位, 屬正向圖法.

///TM2: 二度分帶座標.
///台灣地區中央子午線為東經121度, 投影原點向西平移250,000公尺, 是為Y軸,赤道為X軸.
///澎湖、金門及馬祖等地區, 中央子午線定於東經119度, 投影原點向西平移250000公尺, 是為Y軸, 赤道為X軸.

/// 座標系統          中央經線       東距(X軸)      縮尺系數    扁平率
/// 6度分帶(UTM)   123°(117°)   500,000公尺   0.9996         1/297
/// 3度分帶(TM3)   121°(118°)   350,000公尺   1.0000         1/298.25
/// 2度分帶(TM2)   121°(119°)   250,000公尺   0.9999         1/298.25
/// </summary>
/// <param name="Long"></param>
/// <param name="Lat"></param>
/// <returns>X means TM2 Long, Y means TM2 Lat</returns>
public static DatumPos TWD67EllipsoidaltoTM2(double Long, double Lat)
{
    string[] NewArray = new string[2];
    double UTMNorthing, UTMEasting;
    double N, T, C, A, M;
    double LatRad = ReturnRadian(Lat);
    double LongRad = ReturnRadian(Long);
    double LongOriginRad = ReturnRadian(121); //中央子午線

    //sin^2(x)
    double sin2 = Math.Sin(LatRad) * Math.Sin(LatRad); //**
    //cos^2(x)
    double cos2 = cos2 = Math.Cos(LatRad) * Math.Cos(LatRad); //**
    //q = (a^2 * cos^2(latitude) + b^2 * sin^2(latitude))
    double q = Math.Pow(TWD67a, 2) * cos2 + Math.Pow(TWD67b, 2) * sin2;
    //N = a^2 / q*1/2
    N = Math.Pow(TWD67a, 2) / Math.Sqrt(q);
    //T = tan^2(x) = (1 / cos^2(x)) - 1
    //T = (1/cos2) - 1;
    T = Math.Tan(LatRad) * Math.Tan(LatRad); //**
    //C = (e')^2*cos^2(x)
    C = eccPrimeTWD67Squared * cos2;
    //A = cos(x) * long - 121;
    A = Math.Cos(LatRad) * (LongRad - LongOriginRad);

    M = TWD67a * ((1 - eccTWD67Squared / 4 - 3 * Math.Pow(eccTWD67Squared, 2) / 64 - 5 * Math.Pow(eccTWD67Squared, 3) / 256) * LatRad
        - (3 * eccTWD67Squared / 8 + 3 * Math.Pow(eccTWD67Squared, 2) / 32 + 45 * Math.Pow(eccTWD67Squared, 3) / 1024) * Math.Sin(2 * LatRad)
        + (15 * Math.Pow(eccTWD67Squared, 2) / 256 + 45 * Math.Pow(eccTWD67Squared, 3) / 1024) * Math.Sin(4 * LatRad)
        - (35 * Math.Pow(eccTWD67Squared, 3) / 3072) * Math.Sin(6 * LatRad));

    UTMEasting = OffsetX + (k0 * N * (A + (1 - T + C) * Math.Pow(A, 3) / 6 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeTWD67Squared) * Math.Pow(A, 5) / 120));

    UTMNorthing = (k0 * (M + N * Math.Tan(LatRad) * (Math.Pow(A, 2) / 2 + (5 - T + 9 * C + 4 * C * C) * Math.Pow(A, 4) / 24
        + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeTWD67Squared) * Math.Pow(A, 6) / 720)));

    return new DatumPos()
    {
        X = UTMEasting,
        Y = UTMNorthing
    };
}

 

5. 最後, 再利用 TWD67 二度分帶轉 TWD97 二度分帶座標


private static double TWD67toTWD97_B = 0.000006521;
/// <summary>
/// refer to http://gis.thl.ncku.edu.tw/coordtrans/coordtrans.aspx
/// </summary>
/// <param name="n"></param>
/// <param name="e"></param>
/// <returns></returns>
public static DatumPos TWD67TM2toTWD97TM2(double E67, double N67)
{
    return new DatumPos()
    {
        X = (E67 + 807.8 + TWD67toTWD97_A * E67 + TWD67toTWD97_B * N67),
        Y = (N67 - 248.6 + TWD67toTWD97_A * N67 + TWD67toTWD97_B * E67),
    };
}

 

於是, 我們終於寫出了 WGS84 經緯度轉換成 TWD97 二度分帶的程式:

(後來發現該網頁用的是 TWD67 的二度分帶座標, 所以補上 WGS84 轉 TWD67)


{
    DatumPos pos = WGS84toCartesian(Convert.ToDouble(Long), Convert.ToDouble(Lat), 0);
    pos = SevenParameter_BursaWolf_WGS84_To_TWD67(pos);
    pos = TWD67toEllipsoidal(pos);
    pos = TWD67EllipsoidaltoTM2(pos.X, pos.Y);
    return TWD67TM2toTWD97TM2(pos.X, pos.Y);
}

public static DatumPos WGS84toTWD67TM2(double Long, double Lat)
{
    DatumPos pos = WGS84toCartesian(Convert.ToDouble(Long), Convert.ToDouble(Lat), 0);
    pos = SevenParameter_BursaWolf_WGS84_To_TWD67(pos);
    pos = TWD67toEllipsoidal(pos);
    return TWD67EllipsoidaltoTM2(pos.X, pos.Y);
}

補上 DatumPos 的定義:


/// 表示三度空間的位置
/// </summary>
public class DatumPos
{
    /// <summary>
    /// X , 水平座標值
    /// </summary>
    public double X { get; set; }

    /// <summary>
    /// Y , 垂直座標值
    /// </summary>
    public double Y { get; set; }

    /// <summary>
    /// Z , 高度座標值
    /// </summary>
    public double Z { get; set; }
}

上面的程式可以透過
http://gis.thl.ncku.edu.tw/coordtrans/coordtrans.aspx 驗證正確性,
有一些些誤差就是了.

終於, 我們可以根據位置向台北市政府的網站要一點停車場的資訊了,
從台北市停車管理工程處網頁開始 http://163.29.129.99/pweb/ , 用 IE 搞了半天 (裝一堆外掛啊!)
你應該可以看到

image

點了右上角的 “範圍內停車場查詢”,再點地圖上的某個點, 會看到一個網頁:

image

看到上方的 URL , 我不禁高興了起來, 因為我們已經可以知道 CircleX, CircleY ,
點了距離 (是說, 為什麼預設 Radius = 0 ? 難道這樣查得到停車場嘛? 再點一次距離會讓使用者很高興?)
終於有看到內容…
 

 

 

image

我要做的事情就是把這段 IE 都跑得很辛苦的事情, 交給 Windows Mobile 手機來做.
因為如果你在外面開車, 要找停車位時還要打開 NoteBook, 連上 3.5G 無線網路 (搞不好還是透過手機哩)
然後再打開 IE … 肯定會被旁邊的女性友人笑掉大牙啊!
當然一定是要很帥氣地拿出美美的手機, 隨手點兩下, 停車場資訊就跑到你手指下,
這樣才像是比爾蓋茲叔叔說的 “Information at your fingertips!”呀!

好, 先測一下網址: http://163.29.129.99/?CircleX=304789&CircleY=2772172&Radius=0
image

嗯, 好一個 IIS 建構中畫面 …
這代表寫該網頁的人不想讓我們隨便連上去???
難道使用 Browser , Http 可以騙過工程師嘛…
拿出工具 : Fiddler , 你就會知道, http protocol 根本是光著屁股在街上跑的傢伙,
如果寫網頁的人認為用一些手段 (像是鎖右鍵, 改 Browser Title …等等) 可以讓別人不知道真正的內容,
那還真是國王的新衣耶…

祭出 Fiddler , 重跑台北市停車管理工程處網頁 , 重點的封包內容如下:
image

你會看到, 真正的網址其實是像這樣 http://163.29.129.99/pweb/CirclePK_form.asp?CircleX=308211.868773846&CircleY=2775515.40948218&Radius=1500

然後, 他回傳了 Cookie : ASPSESSIONACTACTAS= …. , 而確切的停車場內容卻是在
image

也就是說, 有了 cookie , 再去跟 http://163.29.129.99/pwb/CirclePK.asp 要內容就有了!
先定義每個停車場的資訊:
 


{
    public string name { get; set; }
    public string space { get; set; }
    public string address { get; set; }
}

取得停車場資訊的程式 (根據實測的結果, 座標應該是 TWD67 二度分帶才是正確的, 所以我做了些修正):


{
    // 將 WGS84 經緯度轉為 TWD67 二度分帶
    DatumPos TM2Pos = Datum.WGS84toTWD67TM2(longitude, latitude);

    CookieContainer cc = new CookieContainer();
    // 取 cookie
    WebUtil.GetPage("http://163.29.129.99/pweb/CirclePK_form.asp?CircleX=" + TM2Pos.X.ToString() +
        "&CircleY=" + TM2Pos.Y.ToString() + "&Radius=" + distanceInMeter.ToString(), null, cc, false, Encoding.GetEncoding("Big5"));

    // 取停車場資訊
    XmlDocument doc = GetXmlDoc(WebUtil.GetPage("http://163.29.129.99/pweb/CirclePK.asp", null, cc, true, Encoding.GetEncoding("Big5")));

    List<ParkInfo> list = new List<ParkInfo>();

    XmlNodeList nodelist = doc.SelectNodes("//tr[@class='bg1']");
    foreach (XmlNode node in nodelist)
    {
        XmlNodeList tdnodes = node.SelectNodes("td");
        if (tdnodes.Count >= 3)
        {
            list.Add(new ParkInfo()
            {
                name = tdnodes[0].InnerText,
                space = tdnodes[1].InnerText,
                address = tdnodes[2].InnerText,
            });
        }
    }

    return list;
}

private static XmlDocument GetXmlDoc(string pageContent)
{
    using (SgmlReader sr = new SgmlReader())
    using (StringReader sReader = new StringReader(pageContent.Trim()))
    {
        try
        {
            sr.InputStream = sReader;
            sr.CaseFolding = Sgml.CaseFolding.ToLower;
            sr.DocType = "HTML";

            XmlDocument xDoc = new XmlDocument();
            xDoc.Load(sr);
            return xDoc;
        }
        catch (XmlException)
        {
            return null;
        }
    }
}

嗯, 你會發現那個處理 HTML 內容的方法, 是透過 SgmlReader 轉為 Xml ,
接下來用 XPath 就可以查到文字內容了.
關於 SgmlReader , 可以參考 微軟釋出 Library : http://code.msdn.microsoft.com/SgmlReader
但是裡面有 bug …後來可以查到又有善心人改的新版 : http://developer.mindtouch.com/SgmlReader
當然, 故事不是這樣就結束了,
我們還必須抓原始碼回來修改符合 .NET CF , 主要修訂內容我就寫一點在這邊,
如果再深入討論 Sgmlreader , 相信會讓大家有點沒力…
這牽扯到 .NET CF 的 Xml Serialization 還是少實做了 SerializationInfo , StreamingContext 的方法.
這兩個是進階的 Xml Serialization 議題, .NET CF 一直到了 3.5 也沒作這個,
所以我們遇到這種情況也只能歎一口氣, 取消這部份 serialization 的功能,
但這點對於 SgmlReader 處理 html 沒有影響, 只是讓部份的 xml serialization /deserialization 遺漏部份資訊而已.

回到主題, 我們開始拖拉 UI , 放一個 ListView, 把 Columns 設兩個 : 停車場, 空位, 然後把寬度也設好,
再把 View 設為 Details, 如下圖:

image

那麼, 取得停車場資訊顯示到 UI 的程式碼如下 (我把距離設定為 1500 公尺, 因為腿短如我, 停車跑不遠…):


{
    GpsPosition pos = gps.GetPosition();
    if ((pos.LatitudeValid == false) || (pos.LongitudeValid == false))
    {
        MessageBox.Show("gps information not ready.");
        return;
    }

    // use a thread to get info!
    Thread t = new Thread(() =>
    {
        var plist = GetParkingInfo(pos.Longitude, pos.Latitude, 1500);
        this.Invoke(new delegateinvoke(() =>
        {
            listView1.Items.Clear();
            foreach (var p in plist)
            {
                listView1.Items.Add(new ListViewItem(new string[] {
                    p.name, p.space}));
            }
        }));
    });
    t.Start();
}

測試運作的結果畫面:

image 

這樣就結束了? 阿…如果你是使用者, 看到這樣的停車場資訊, 會嫌不夠吧!
所以當然還要補上, 當使用者點選其中的停車場時,
顯示更多的資訊啊!
我們可以透過類似這樣的 URL 取得資訊 http://163.29.129.99/pweb/QueryPkDetail.asp?PKName=%C0%F4%A4s%A5%AD%AD%B1%B0%B1%A8%AE%B3%F5

於是我們再新增一個 Form, 放個 label,
在 Listview 的 ItemActivate event 寫下程式碼:


{
    if (listView1.SelectedIndices.Count > 0)
    {
        ListViewItem firstli = listView1.Items[listView1.SelectedIndices[0]];

        string pkname = firstli.Text;

        // show wait cursor
        Cursor.Current = Cursors.WaitCursor;

        Encoding big5 =  Encoding.GetEncoding("Big5");

        XmlDocument doc = GetXmlDoc(WebUtil.GetPage("http://163.29.129.99/pweb/QueryPkDetail.asp?PKName="+
            HttpUtility.UrlEncode(pkname, big5), null, null, true, big5));

        StringBuilder displaycontext = new StringBuilder();
        displaycontext.AppendLine(pkname);

        XmlNodeList contentnodes = doc.SelectNodes("//td[@class='ExMenu']");
        
        for (int i=0;i<contentnodes.Count-1;i++)
        {
            if (i == 0)
                displaycontext.AppendLine("空位資訊:" + contentnodes[0].InnerText);
            else
                displaycontext.AppendLine(contentnodes[i].InnerText);
        }

        // reset cursor to default.
        Cursor.Current = Cursors.Default;

        using (ParkingDetailForm pkd = new ParkingDetailForm())
        {
            pkd.DisplayInfo(displaycontext.ToString());
        }

    }
}

你應該可以注意到,
這次利用 Cursor.Current = Cursors.WaitCursor; 來顯示漏斗游標,
告訴使用者說, 程式正在努力下載停車場詳細資訊.

執行結果如下:
image

成功了, We did it!  登登登登~~~ (Dora 看太多了啊~~~)

原始碼下載  : TaipeiParking.zip