0

我已经用网格点建模了一个3D 地球，如下所示：

地球

这些点在 3D 空间中表示为XYZ坐标。然后我根据我从这里获取的脚本 将XYZ转换为Lat/Long/Height （海拔）： JSFiddle

出于某种原因，当我试图从我的集合中查找LLH的XYZ时，我得到了非常奇怪的结果，因此我尝试通过将XYZ转换为LLH，然后将相同的LLH 转换回XYZ来验证初始脚本，以查看我是否得到相同的坐标。

相反，结果坐标是地球上的某个XYZ ，与原始XYZ位置无关。

XYZ 到 LLH 脚本：
来源：JSFiddle

def xyzllh(x,y,z):
    """        xyz vector  to  lat,lon,height
         output:
    
              llhvec[3] with components
    
              flat      geodetic latitude in deg
              flon      longitude in deg
              altkm     altitude in km
    """
    dtr =  math.pi/180.0
    rrnrm = [0.0] * 3
    llhvec = [0.0] * 3
    
    geodGBL()
    
    esq = EARTH_Esq

    rp = math.sqrt( x*x + y*y + z*z )
    flatgc = math.asin( z / rp )/dtr
    testval= abs(x) + abs(y)
    
    if ( testval < 1.0e-10):
        flon = 0.0
    else:
        flon = math.atan2( y,x )/dtr 
    if (flon < 0.0 ):
        flon = flon + 360.0
    
    p = math.sqrt( x*x + y*y )
    
    # on pole special case
    
    if ( p < 1.0e-10 ):
        flat = 90.0
        if ( z < 0.0 ):
            flat = -90.0
        altkm = rp - rearth(flat)
        llhvec[0]  = flat
        llhvec[1]  = flon
        llhvec[2]  = altkm
        return  llhvec

    # first iteration, use flatgc to get altitude 
    # and alt needed to convert gc to gd lat.
    
    rnow = rearth(flatgc)
    altkm = rp - rnow
    flat = gc2gd(flatgc,altkm)
    
    rrnrm = radcur(flat)
    rn = rrnrm[1]
    
    for x in range(5):
        slat = math.sin(dtr*flat)
        tangd = ( z + rn*esq*slat ) / p
        flatn = math.atan(tangd)/dtr
    
        dlat = flatn - flat
        flat = flatn
        clat = math.cos( dtr*flat )
    
        rrnrm = radcur(flat)
        rn = rrnrm[1]
    
        altkm = (p/clat) - rn
    
        if ( abs(dlat) < 1.0e-12 ):
            break

    llhvec[0]  = flat
    llhvec[1]  = flon
    llhvec[2]  = altkm
    return llhvec 
# globals
EARTH_A = 0
EARTH_B = 0
EARTH_F = 0
EARTH_Ecc = 0
EARTH_Esq = 0

# starting function do_llhxyz()
CallCount = 0
llh  = [0.0] * 3
dtr = math.pi/180

CallCount = CallCount + 1

sans = " \n"

llh = xyzllh(x,y,z)

latitude = llh[0]
longitude= llh[1]
hkm = llh[2]
height = 1000.0 * hkm

latitude = fformat(latitude,5)
longitude = fformat(longitude,5)
height = fformat(height,1)

sans = sans +"Latitude,Longitude, Height (ellipsoidal) from ECEF\n"
sans = sans + "\n"
sans = sans +"Latitude  : " + str(latitude)  + "   deg N\n"
sans = sans +"Longitude : " + str(longitude - 180) + "   deg E\n"
sans = sans +"Height    : " + str(height)    + "   m\n"

lats = []
longs = []
heights = []
lats.append(str(latitude))
longs.append(str(longitude - 180))
heights.append(str(height))

这是 LLH 到 XYZ 的脚本：
来源：www.mathworks.com

a = 6378137
t = 8.1819190842622e-2

# (prime vertical radius of curvature)
N = a / math.sqrt(1 - (t*t) * (math.sin(lat)*math.sin(lat)))

x = []
y = []
z = []

# results:
x.append( ((N+height) * math.cos(lat) * math.cos(long))/1000 )
y.append( ((N+height) * math.cos(lat) * math.sin(long))/1000 )
z.append( (((1-t*t) * N + height) * math.sin(lat))/1000 )

有人知道我在这里做错了什么吗？谢谢！

4

0 回答 0