UTM to lt/lg converting problem

Help
dou wen
2005-09-16
2013-04-17
  • dou wen
    dou wen
    2005-09-16

    I have tested the converting from UTM coordinate to normal longitude/latitude of qpegps,
    the converting stuff is mainly based on following two methods of Class MapUTM, and i found something wrong with it.

    for example, use following qpegps code, the UTM  coordinate (651731E,3266424N,48R) will be converted to
    (29.52510N,105.78291E), not the correct answer (29.51822N,106.56559E), i am not familiar with algorithm detail for this converting, but i think this is a little bug of qpegps.

    the test is based on qpegps 0.9.2.3 pre

    thanks.

    --------------------------------------------------------------------

    MapUTM::MapUTM(QString *mapInfo, bool utmCoord):MapBase(mapInfo)
    {
    ...
    ...
        xlong1 = 0.5*atanh(cos(latitude1) * sin(checkLg(longitude1-stdLong)));
        ylat1 = atan(tan(latitude1) / cos(checkLg(longitude1-stdLong)));
        xlong2 = 0.5*atanh(cos(latitude2) * sin(checkLg(longitude2-stdLong)));
        ylat2 = atan(tan(latitude2) / cos(checkLg(longitude2-stdLong)));
    ...
    ...
    }

    bool MapUTM::calcltlg(double *lt, double *lg, double x, double y)
    {
        double xt,yt;
        xt = xlong1 + ( (x - x1) * checkLg(xlong2 - xlong1) / (x2-x1) );
        yt = ylat1 + ( (y - y1) * (ylat2 - ylat1) / (y2-y1));
                                                                                   
        printf("orig xlong1 = %f\n",xlong1);
        printf("orig ylat1 = %f\n",ylat1);
                                                                                   
        *lt = asin(sin(yt)/cosh(xt));
        *lg = checkLg(stdLong + atan(sinh(xt) / cos(yt)));
        return TRUE;
    }

     
    • The mentioned 2 methods do the conversion between longitude/latitude and the corresponding pixels on the map for a transverse mercator projection (UTM is a transverse mercator projection where the coordinates are giving in northing,easting and zone instead of longitude and latitude). The transformation from UTM coordinates to latitude/longitude is done with the method UTMtoLL. As I do not have a single UTM map, I couldn't check it myself if this works correct. Another important point of this conversion, is on which "datum" your map is based.
      Are the UTM coords and your correct answer both based on the WGS84 datum ? The current algorithm works only for this datum. Is there a common need for UTM maps with another datum than WGS84?

       
    • dou wen
      dou wen
      2005-09-19

      hi, Ralf, of course you are right ,the UTMtoLL method is used  when importing a new UTM map, the MapUTM constructor will call 
      the UTMtoLL to get the longitude/latitude from input 2 reference point whose  coordinates are giving in northing,easting and zone
      .  i am sure and checked the caculation in UTMtoLL is all right. so i didnt mention this method.

      _but_   the MapUTM constructor will also pre-caculate xlong1, xlat1,xlong2,xlat2, which will be used in the MapUTM::calcltlg method to caculate  the lt/lg value when i  mouse-click any point of the map.

      i.e.,  when i mouse click any point on this imported map
      , the caculating for the lt/lg coordinates (in calcltlg) is not depend on
      longitude1/latitude1 and  longitude1/latitude2 of the two reference point  directly, but it depend directly on xlong1/xlat1 and xlong2/xlat2 which were "pre-calculated " when the map opject was constructed. 
        I think the caculation for the lt/lg coordinates is not correct.

        That is the reason i list only the MapUTM constructor and calcltlg methods.
      here is my guess:
      (1) i guess the pre-calc for xlong1/xlat1 and xlong2/xlat2 is not correct.
      (2) i guess the algorithm in MapUTM::calcltlg is not correct.

      I am very sure the map i mentioned is WGS84 datum based UTM map.

      Ralf,  here is a UTM map item of my maps.txt, you could simply generate a same size png file of the map item and import it to
      your test system to  verify the problem what i prompt.

      "UTM cq.png 10000 4134 5984 48R 3279760 639498 236 72 3266352 650470 2900 3581"

      you can also use a public UTM2lat converter to verify what i said ( http://sinfogeo.com/tools/utm2lat.htm  )  .

      i would also like to send the OZI explorer map to you if you accept.

      Best regards

       
      • Hi Dou,

        sorry I wasn't reading your post carefully enough.
        Thanks for the UTM string (wow, this map is really big; please send it to rhaselme@users.sourceforge.net).

        The projection algorithm works this way:
        UTMtoLL: The UTM coordinates are converted to lat/long
        calcxy: converts lat/long to pixel coordinates (=x,y)
        calcltlg: conversts pixels to lat/long
        The projection calculates coordinates in the range of -1 to +1( resp from-PI to +PI ..., depends on projection). So there is always a linear interpolation from this (let's call it "intermediate linear coord system"=ilcs) to the coordinates in pixels.
        Precalc in constructor: calculates the values of the 2 reference points in the "ilcs"(xlong1/2,ylat1/2).

        The projection equations are:
        x = 0.5 * atanh(B)
        y = atan( tan(lat) / cos(long -stdLong) -y0
        inverse:
        lat = asin(sin(D)/cos(x))
        long = stdLong + atan(sinh(x)/cos(D))
        where
        B = cos(lat)*sin(long-stdLong)
        D = y +y0
        As y0 is a common offset, it shouldn't have any impact after the linear interpolation, so I set it to 0.

        The only posible problem I saw is the usage of the checkLg function. So please try this (modified) methods:
        bool MapUTM::calcxy(double *x, double *y, double lg, double lt)
        {
            double xt,yt;
            xt = 0.5*atanh(cos(lt) * sin(checkLg(lg-stdLong)));
            yt = atan(tan(lt) / cos(checkLg(lg-stdLong)));
            *x = x1 + ( xt - xlong1 * (x2-x1) / xlong2 - xlong1);
            *y = y1 + ( (yt - ylat1) * (y2-y1) / (ylat2 - ylat1) );
            if ((*x < 0)||(*x >= mapSizeX)||(*y < 0)||(*y >= mapSizeY))
                return (FALSE);
            return (TRUE);
        }

        bool MapUTM::calcltlg(double *lt, double *lg, double x, double y)
        {
            double xt,yt;
            xt = xlong1 + ( (x - x1) * xlong2 - xlong1 / (x2-x1) );
            yt = ylat1 + ( (y - y1) * (ylat2 - ylat1) / (y2-y1));
            *lt = asin(sin(yt)/cosh(xt));
            *lg = checkLg(stdLong + atan(sinh(xt) / cos(yt)));
            return TRUE;
        }

        Regards,
        Ralf

         
    • dou wen
      dou wen
      2005-09-22

      hi,Ralf
      i have tried the two modified methods, but the test result is still
      not correct :-(

      Best regards