自从接触mapinfo,被它的易用性打动,可惜国外产的软件,没有咱们国家的等差分纬线多圆锥投影,所以很多专题地图(世界的)要在mapinfo中描图的话就不能与其他数据图层叠加,这个投影找了好久,终于今天看到了它的正反解源程序,相信这个问题很快会得到解决,如果成功的话,那mapinfo就更完美了,更适合我们制作空白图,为了这个目标奋斗,先兴奋一下,其实问题也不少。

'正切差分纬线多圆锥投影的正解变换public sub latlontoxysecant(lat as double, lon as double, x as double, y as double)    dim latt as double, lont as double    dim y0 as double, r as double, u0 as double    dim yn as double, xn as double    dim deltaqn as double, deltaqi as double    dim lonn as double, u0r100 as double, rou as double    dim yn_y0 as double

    latt = abs(lat) * cd    lont = lon * cd

    r = 6371116    u0 = 1# / 14000000    u0r100 = 45.5079714285714 ' u0 * r * 100    lonn = 210 * cd        y0 = (latt 0.06683225 * (latt ^ 4)) * u0r100    yn = y0 0.20984 * latt * u0r100    yn_y0 = 0.20984 * latt * u0r100    if (abs(yn) > 112) then        yn = 112    end if    xn = sqr(112 * 112 - yn * yn) 20    if (lat = 0) then '赤道        y = 0        x = xn / lonn * (1.1 - 0.11106126 * tan(lont / 5)) * lont    else        rou = (xn ^ 2 yn_y0 ^ 2) / 2 / yn_y0        deltaqn = arcsin(xn / rou)            deltaqi = deltaqn / lonn * (1.1 - 0.11106126 * tan(lont / 5)) * lont        y = y0 rou * (1 - cos(deltaqi))        x = rou * sin(deltaqi)    end if    if (lat < 0) then        y = -y    end ifend sub

'正切差分纬线多圆锥投影反解变换方法public sub xytolatlonsecant(xc as double, yc as double, lat as double, lon as double)    dim r as double, u0 as double    dim y0 as double, y01 as double    dim xn as double, yn as double    dim xn1 as double, yn1 as double    dim rou as double, rou1 as double    dim yn_y0 as double    dim u0r100 as double, n as integer    dim lat1 as double, lon1 as double, lonn as double    dim deltaqi as double, deltaqi1 as double, deltaqn as double    dim x as double, y as double, yy as double    dim f as double, f1 as double, f2 as double

    x = (xc - wcsx0) / xyfact    y = (wcsy0 - yc) / xyfact

    on error resume next

    r = 6371116    u0 = 1# / 14000000    lonn = 210 * cd    u0r100 = 45.5079714285714 ' u0 * r * 100        '反解纬度    lat = abs(y / u0r100)        if (lat > pi / 4) then        lat = pi / 4    end if    n = 0    do        n = n 1

        y0 = (lat 0.06683225 * (lat ^ 4)) * u0r100        yn = y0 0.20984 * lat * u0r100        xn = sqr(112 * 112 - yn * yn) 20                    yn_y0 = 0.20984 * lat * u0r100                rou = (xn * xn yn_y0 ^ 2) / 2 / yn_y0        y01 = (1 0.267329 * (lat ^ 3)) * u0r100                    yn1 = y01 9.5493        xn1 = -yn * yn1 / sqr(112 * 112 - yn * yn)        rou1 = (2 * yn_y0 * xn * xn1 - xn * xn * 9.5493) / yn_y0 ^ 2 / 2 9.5493 / 2

        f = (y0 rou - abs(y)) ^ 2 x * x - rou * rou        f1 = 2 * (y0 rou - abs(y)) * (y01 rou1) - 2 * rou * rou1        f2 = f / f1                        lat1 = lat - f2        if (abs(lat1 - lat) < 0.00000001 or n > 100 or abs(f) < 0.00000001) then exit do        lat = lat1    loop

    '反解经度    lon = x / u0r100    n = 0    if (y = 0) then '在赤道上        xn = 112 20        do            n = n 1                        f = xn / lonn * (1.1 - 0.11106126 * tan(lon / 5)) * lon - x            f1 = xn / lonn * (1.1 - 0.11106126 * (lon / 5 / (cos(lon / 5) ^ 2) tan(lon / 5)))            f2 = f / f1            lon1 = lon - f2            if (abs(lon1 - lon) < 0.00000001 or n > 100 or abs(f) < 0.00000001) then                exit do            end if            lon = lon1        loop    else        do            n = n 1            y0 = (lat 0.06683225 * (lat ^ 4)) * u0r100            yn = y0 0.20984 * lat * u0r100            yn_y0 = 0.20984 * lat * u0r100                        xn = sqr(112 * 112 - yn * yn) 20

            rou = (xn * xn yn_y0 ^ 2) / 2 / yn_y0            deltaqn = arcsin(xn / rou)            deltaqi = arcsin(x / rou)                        f = deltaqi - deltaqn / lonn * (1.1 - 0.11106126 * tan(lon / 5)) * lon            f1 = -deltaqn / lonn * (1.1 - 0.11106126 * tan(lon / 5) - 0.02221225 * lon / (cos(lon / 5) ^ 2))            f2 = f / f1            lon1 = lon - f2            if (abs(lon1 - lon) < 0.00000001 or n > 100 or abs(f) < 0.00000001) then                exit do            end if            lon = lon1        loop    end if    if (y < 0) then        lat = -lat    end if    lon = lon / cd    lat = lat / cdend sub