簡介坐標(biāo)轉(zhuǎn)換是測繪與地理信息領(lǐng)域高頻且基礎(chǔ)的技術(shù)環(huán)節(jié)尤其是存量西安80數(shù)據(jù)與新一代CGCS2000成果并存時(shí)可靠轉(zhuǎn)換工具尤為關(guān)鍵。資源包正是一套聚焦CGCS2000與西安80互轉(zhuǎn)的精準(zhǔn)轉(zhuǎn)換工具包面向測繪、地理信息及工程測量人員解決新舊坐標(biāo)系統(tǒng)間的數(shù)據(jù)兼容與精度統(tǒng)一問題。壓縮包共97個(gè)文件合計(jì)約4.26MB核心含COORD GM2.0轉(zhuǎn)換程序exe、49個(gè)坐標(biāo)參數(shù)文件cod、CSV樣例數(shù)據(jù)、TXT參數(shù)說明、ZGF/GGF格網(wǎng)文件及GIF/JPG操作截圖支持七參數(shù)、三參數(shù)等常用轉(zhuǎn)換與批量處理。目前已有3152人學(xué)習(xí)下載尤其適合需要坐標(biāo)轉(zhuǎn)換、參數(shù)校核以及系統(tǒng)學(xué)習(xí)坐標(biāo)系統(tǒng)間換算流程的從業(yè)者和相關(guān)專業(yè)學(xué)生。通過本包既能直接運(yùn)行軟件完成坐標(biāo)互轉(zhuǎn)也能結(jié)合參數(shù)文件理解橢球、基準(zhǔn)面與投影變換快速應(yīng)用到實(shí)際測繪或GIS項(xiàng)目中。1. 國家大地2000精準(zhǔn)轉(zhuǎn)換器坐標(biāo)轉(zhuǎn)換不是換個(gè)數(shù)字是一套參數(shù)和模型的事只要做過一輪測繪數(shù)據(jù)匯交基本都遇到過“坐標(biāo)轉(zhuǎn)換”這四個(gè)字。國家大地2000CGCS2000坐標(biāo)系已經(jīng)是當(dāng)前成果入庫和國土空間規(guī)劃數(shù)據(jù)庫的默認(rèn)基準(zhǔn)手里那批北京54、西安80、WGS84甚至地方獨(dú)立坐標(biāo)的存量數(shù)據(jù)必須經(jīng)過“國家大地2000精準(zhǔn)轉(zhuǎn)換器”才能并網(wǎng)。很多項(xiàng)目卡在這一步是因?yàn)檗D(zhuǎn)換不只是把度分秒換成米而是橢球、基準(zhǔn)、投影、七參數(shù)一路都要對。這篇內(nèi)容就按一線數(shù)據(jù)處理的順序把坐標(biāo)系理論、pyproj轉(zhuǎn)換器搭建、七參數(shù)擬合和常見返工原因講透。適合正在接數(shù)據(jù)匯交任務(wù)的測繪技術(shù)員、GIS開發(fā)者和被存量數(shù)據(jù)折騰的從業(yè)者。2. 轉(zhuǎn)換器底層邏輯橢球、基準(zhǔn)、投影三個(gè)東西先分清2.1 橢球、基準(zhǔn)、投影三個(gè)概念決定轉(zhuǎn)換器怎么寫做坐標(biāo)轉(zhuǎn)換最容易出問題的地方就是把“橢球”“基準(zhǔn)”“投影”混為一談。橢球是一個(gè)數(shù)學(xué)曲面北京54用克拉索夫斯基橢球西安80用IAG75橢球CGCS2000用2000國家大地橢球橢球定了經(jīng)緯度才有依附的對象。基準(zhǔn)是橢球與地球?qū)嶋H位置之間的固定關(guān)系它決定某個(gè)點(diǎn)的經(jīng)緯度對應(yīng)地球上的哪塊地。投影則是把球面坐標(biāo)展開成平面坐標(biāo)的算法。轉(zhuǎn)換器正是在這三個(gè)層面工作同一基準(zhǔn)下可以做“經(jīng)緯度換平面坐標(biāo)”的投影計(jì)算不同基準(zhǔn)之間則需要先做基準(zhǔn)轉(zhuǎn)換最常見的就是七參數(shù)或四參數(shù)。我常見到有人拿CGCS2000的投影參數(shù)去處理西安80的平面坐標(biāo)出來的結(jié)果偏移幾百米原因就是跳過了基準(zhǔn)轉(zhuǎn)換這一層。所以寫轉(zhuǎn)換器之前先回答三個(gè)問題源數(shù)據(jù)是哪種坐標(biāo)系它屬于哪個(gè)投影帶目標(biāo)坐標(biāo)到底是經(jīng)緯度還是高斯平面坐標(biāo)這三個(gè)問題答不清楚后面所有代碼都是白寫。2.2 CGCS2000 與 WGS84、北京54、西安80 的差異CGCS2000是國家大地坐標(biāo)系的原點(diǎn)位于地球質(zhì)心與ITRF國際地球參考框架在2000歷元一致。它的橢球長半軸a6378137米扁率f1/298.257222101。單看橢球參數(shù)WGS84的扁率是1/298.257223563兩者在地表引起的橢球面差不足0.1毫米這也是很多人問“WGS84能不能直接當(dāng)CGCS2000用”的原因。但實(shí)際不能直接畫等號。兩套坐標(biāo)系雖然橢球幾乎一樣框架對準(zhǔn)方式和觀測歷元并不完全相同不同時(shí)期的WGS84數(shù)據(jù)相對CGCS2000可能存在幾厘米到幾分米的系統(tǒng)性偏差。對于毫米級控制網(wǎng)這個(gè)偏差不能忽略對于普通不動產(chǎn)權(quán)籍?dāng)?shù)據(jù)有的地方直接沿用也能通過檢查但我一般不建議這么干。正確做法是先拿項(xiàng)目區(qū)域內(nèi)的已知點(diǎn)做比對再決定是否需要參數(shù)轉(zhuǎn)換。北京54和西安80則完全不同。它們屬于參心坐標(biāo)系橢球中心不與地心重合與CGCS2000之間是整體平移加旋轉(zhuǎn)的關(guān)系。這個(gè)偏移通常在幾十米到上百米必須用七參數(shù)或四參數(shù)做基準(zhǔn)轉(zhuǎn)換不能靠加減固定數(shù)解決。實(shí)際操作中每個(gè)區(qū)域的參數(shù)都不一樣這就是為什么不能隨便下載一個(gè)“全國通用參數(shù)”到處套。2.3 高斯-克呂格投影與帶號換算拿到米制坐標(biāo)前先做的一道題CGCS2000平面坐標(biāo)普遍采用高斯-克呂格投影分3度帶和6度帶兩種。6度帶從中央子午線135°E逆推帶號與中央子午線的關(guān)系是L6n-33度帶則是L3n。我國大陸主要使用3度帶尤其是1:1萬及更大比例尺數(shù)據(jù)。中央子午線定錯(cuò)會讓平面坐標(biāo)的Y方向出現(xiàn)幾百公里級的偏移這種錯(cuò)誤在數(shù)據(jù)檢查中最常見。以3度帶為例37帶的中央子午線是111°E38帶是114°E39帶是117°E。一個(gè)點(diǎn)落在經(jīng)度116.4°E附近按常規(guī)應(yīng)使用39帶117°E中央子午線或自定義中央經(jīng)線而不應(yīng)該套用38帶。很多老數(shù)據(jù)的帶號信息丟了轉(zhuǎn)換前只能根據(jù)坐標(biāo)Y值的數(shù)量級反推——Y坐標(biāo)大約在500公里上下浮動如果看到Y(jié)是38開頭或39開頭基本能判斷帶號。投影參數(shù)中還有一個(gè)容易被忽略的點(diǎn)x_0東偏移量通常設(shè)為500000米是為了避免中央子午線以西出現(xiàn)負(fù)坐標(biāo)。有些地方獨(dú)立坐標(biāo)系會把x_0改成自定義值轉(zhuǎn)換前必須核對原始坐標(biāo)系定義否則成果直接不能用。3. 用 pyproj 搭建 CGCS2000 轉(zhuǎn)換器最小可運(yùn)行代碼與參數(shù)解讀3.1 搭建環(huán)境并完成首次CGCS2000投影轉(zhuǎn)換做坐標(biāo)轉(zhuǎn)換我一般直接用pyproj庫它是PROJ的Python綁定EPSG代碼庫和坐標(biāo)變換模型都內(nèi)置比手寫投影公式可靠得多。安裝命令非常簡單pip install pyproj裝好后下面這段代碼把WGS84經(jīng)緯度轉(zhuǎn)成CGCS2000高斯平面坐標(biāo)。這里用的是3度帶中央子午線117°E的參數(shù)from pyproj import CRS, Transformer # CGCS2000經(jīng)緯度坐標(biāo)EPSG:4490 crs_geo_2000 CRS.from_epsg(4490) # 中央子午線117°E的CGCS2000高斯-克呂格投影 crs_tm_117 CRS.from_proj4( projtmerc lat_00 lon_0117 k1 x_0500000 ellpsGRS80 unitsm no_defs ) trans Transformer.from_crs(crs_geo_2000, crs_tm_117, always_xyTrue) # 注意always_xyTrue 時(shí)先傳經(jīng)度再傳緯度 lon, lat 116.3972, 39.9075 x, y trans.transform(lon, lat) print(fX{x:.3f}, Y{y:.3f})這段代碼的語義是先定義CGCS2000經(jīng)緯度坐標(biāo)系再定義以117°E為中央子午線的高斯克呂格投影坐標(biāo)系然后通過Transformer完成投影轉(zhuǎn)換。CRS.from_proj4里每個(gè)參數(shù)都有明確含義projtmerc表示高斯-克呂格投影l(fā)at_0是原點(diǎn)緯度lon_0是中央子午線k是比例因子x_0是東偏移量ellpsGRS80對應(yīng)CGCS2000橢球。實(shí)際跑這段代碼前要再次確認(rèn)輸入坐標(biāo)的順序。always_xyTrue表示輸入輸出統(tǒng)一按經(jīng)度在前、緯度在后的順序如果你源數(shù)據(jù)是緯度在前需要先調(diào)整否則結(jié)果會錯(cuò)到另一個(gè)半球。3.2 批量轉(zhuǎn)換CSV坐標(biāo)數(shù)據(jù)從經(jīng)緯度到高斯平面坐標(biāo)生產(chǎn)項(xiàng)目里很少一次只轉(zhuǎn)一個(gè)點(diǎn)更多是幾千上萬條記錄。批量轉(zhuǎn)換時(shí)用pandas讀CSV再用numpy數(shù)組一次性交給Transformer處理效果遠(yuǎn)好于循環(huán)調(diào)用import pandas as pd import numpy as np from pyproj import CRS, Transformer df pd.read_csv(source_points.csv) trans Transformer.from_crs( CRS.from_epsg(4326), # WGS84經(jīng)緯度 CRS.from_proj4( projtmerc lat_00 lon_0120 k1 x_0500000 ellpsGRS80 unitsm no_defs ), always_xyTrue, ) lon df[L].to_numpy() lat df[B].to_numpy() x, y trans.transform(lon, lat) df[X] np.round(x, 3) df[Y] np.round(y, 3) df.to_csv(converted_2000.csv, indexFalse)這里把WGS84經(jīng)緯度直接轉(zhuǎn)到CGCS2000平面坐標(biāo)中間沒有做顯式的基準(zhǔn)轉(zhuǎn)換相當(dāng)于默認(rèn)兩套框架差異在本項(xiàng)目精度要求下可接受。如果你的源數(shù)據(jù)是老GPS設(shè)備采集或者項(xiàng)目要求厘米級精度建議先按第4章的做法用公共控制點(diǎn)校準(zhǔn)一下而不是直接信任這一步。批量轉(zhuǎn)換時(shí)務(wù)必保證源CSV中沒有缺失值和超范圍數(shù)據(jù)。Transformer對空值會直接拋異常生產(chǎn)腳本里最好在前面加一行df df.dropna(subset[L, B])避免跑一半中斷。3.3 中央子午線和帶號的選擇邏輯很多人在這一步翻車是因?yàn)椴磺宄撚媚膫€(gè)中央子午線。我一般會先看源數(shù)據(jù)元數(shù)據(jù)沒有就按坐標(biāo)Y值反推Y在500公里上下說明已經(jīng)處于投影帶內(nèi)結(jié)合X方向量級判斷大概區(qū)域再匹配對應(yīng)的3度帶帶號。還有一種常見情況源數(shù)據(jù)是6度帶目標(biāo)要求3度帶這不能直接帶號相減。正確做法是先投影反算經(jīng)緯度再重新投影到新的3度帶中央子午線。pyproj里可以用兩次Transformer完成或者用PROJ的pipeline功能一步到位。建議第一次做這種跨帶轉(zhuǎn)換時(shí)先在已知點(diǎn)上驗(yàn)證不要一萬條數(shù)據(jù)直接跑完再檢查。4. 精準(zhǔn)的前提是七參數(shù)用公共控制點(diǎn)擬合可靠轉(zhuǎn)換參數(shù)4.1 收集公共控制點(diǎn)與坐標(biāo)制式預(yù)檢坐標(biāo)轉(zhuǎn)換要做到“精準(zhǔn)”最核心的一步是參數(shù)本身可靠。所謂公共控制點(diǎn)就是同一個(gè)物理點(diǎn)在源坐標(biāo)系和目標(biāo)坐標(biāo)系里都有坐標(biāo)的那些點(diǎn)。常見來源是控制網(wǎng)成果表、CORS站已知點(diǎn)、已有竣工測量控制點(diǎn)。我一般要求至少6個(gè)公共點(diǎn)點(diǎn)位盡量覆蓋項(xiàng)目區(qū)四周和中部不要全部集中在一條線上。拿到公共點(diǎn)后先做一次坐標(biāo)制式預(yù)檢哪些點(diǎn)是經(jīng)緯度哪些是平面坐標(biāo)平面坐標(biāo)帶號是多少高程是正常高還是大地高。七參數(shù)計(jì)算需要的是地理坐標(biāo)對應(yīng)的地心直角坐標(biāo)XYZ如果你的控制點(diǎn)是平面坐標(biāo)必須先反算成經(jīng)緯度再按橢球參數(shù)轉(zhuǎn)成XYZ。這一步最容易出問題的是高程體系不統(tǒng)一。坐標(biāo)轉(zhuǎn)換嚴(yán)格說是在三維空間進(jìn)行的如果有的點(diǎn)帶高程、有的點(diǎn)高程缺失后續(xù)殘差會很大。我的習(xí)慣是統(tǒng)一用大地高沒有大地高的點(diǎn)直接剔除。4.2 最小二乘求解布爾沙七參數(shù)七參數(shù)模型里有3個(gè)平移參數(shù)、3個(gè)旋轉(zhuǎn)參數(shù)和1個(gè)尺度參數(shù)。至少需要3個(gè)公共點(diǎn)才能求出唯一解但實(shí)際中公共點(diǎn)多于3個(gè)時(shí)可以用最小二乘平差提高可靠性。下面是求解核心代碼import numpy as np def build_A_l(src_xyz, tgt_xyz): 構(gòu)建布爾沙模型誤差方程矩陣。 src_xyz: 源坐標(biāo)系地心直角坐標(biāo), 形狀(n, 3) tgt_xyz: 目標(biāo)坐標(biāo)系地心直角坐標(biāo), 形狀(n, 3) 返回 A 和 l用于求解 [tx, ty, tz, rx, ry, rz, m] A_rows, l_rows [], [] for (x1, y1, z1), (x2, y2, z2) in zip(src_xyz, tgt_xyz): # X方向方程: x2 - x1 tx z1*ry - y1*rz x1*m A_rows.append([1, 0, 0, 0, z1, -y1, x1]) l_rows.append([x2 - x1]) # Y方向方程: y2 - y1 ty - z1*rx x1*rz y1*m A_rows.append([0, 1, 0, -z1, 0, x1, y1]) l_rows.append([y2 - y1]) # Z方向方程: z2 - z1 tz y1*rx - x1*ry z1*m A_rows.append([0, 0, 1, y1, -x1, 0, z1]) l_rows.append([z2 - z1]) return np.array(A_rows), np.array(l_rows) # 公共點(diǎn)XYZ坐標(biāo) src np.array([ [1_000_000, 2_000_000, 3_000_000], [1_100_000, 2_100_000, 3_100_000], ... ]) tgt np.array([ [1_000_010, 2_000_020, 3_000_030], ... ]) A, l build_A_l(src, tgt) x np.linalg.inv(A.T A) A.T l tx, ty, tz, rx, ry, rz, m x.flatten() print(f平移: {tx:.3f}, {ty:.3f}, {tz:.3f} m) print(f旋轉(zhuǎn): {rx:.6f}, {ry:.6f}, {rz:.6f} 弧度) print(f尺度: {m:.9f})代碼邏輯是先把每個(gè)公共點(diǎn)的三維坐標(biāo)代入布爾沙模型形成誤差方程再把所有方程組裝成矩陣用最小二乘法解出7個(gè)參數(shù)。公式里旋轉(zhuǎn)參數(shù)的單位是弧度尺度參數(shù)m是無量綱的比值。實(shí)際使用時(shí)要特別注意不同文檔里的布爾沙模型旋轉(zhuǎn)方向符號可能相反。拿到別人給的七參數(shù)后我的習(xí)慣是先在一個(gè)已知點(diǎn)上試算確認(rèn)殘差在毫米級再批量用。4.3 殘差檢核與四參數(shù)選擇邊界七參數(shù)求出來后必須把所有公共點(diǎn)代回去算殘差。殘差是目標(biāo)坐標(biāo)與轉(zhuǎn)換結(jié)果之間的差一般按點(diǎn)位中誤差統(tǒng)計(jì)def apply_transform(xyz, params): tx, ty, tz, rx, ry, rz, m params x, y, z xyz.T x2 x tx m * x - rz * y ry * z y2 y ty rz * x m * y - rx * z z2 z tz - ry * x rx * y m * z return np.column_stack([x2, y2, z2]) pred apply_transform(src, x) rmse np.sqrt(np.mean(np.sum((pred - tgt) ** 2, axis1))) print(fRMS: {rmse:.3f} m)殘差如果超過2厘米對控制網(wǎng)級別的轉(zhuǎn)換就不能接受需要回頭檢查公共點(diǎn)坐標(biāo)是否有粗差。實(shí)踐中我通常會留出三分之一公共點(diǎn)不參與求參專門做驗(yàn)證避免把所有點(diǎn)都用于擬合導(dǎo)致“看起來精度很高、實(shí)際一測就偏”。另外如果你的項(xiàng)目區(qū)域很小比如一個(gè)小區(qū)或一片廠區(qū)用七參數(shù)反而容易因?yàn)辄c(diǎn)位分布太集中導(dǎo)致旋轉(zhuǎn)參數(shù)不穩(wěn)定。這種情況下用四參數(shù)更合適四參數(shù)在平面上做平移、旋轉(zhuǎn)、縮放只需要2個(gè)公共點(diǎn)就能解算適合1平方公里內(nèi)的局部轉(zhuǎn)換。5. 坐標(biāo)轉(zhuǎn)換避坑指南五個(gè)最容易讓精度翻車的細(xì)節(jié)5.1 帶號搞錯(cuò)Y坐標(biāo)憑空多出幾百公里現(xiàn)象轉(zhuǎn)換后的坐標(biāo)X還在正常范圍內(nèi)但Y值出現(xiàn)38400000這類數(shù)字一看就是6度帶和3度帶混用了。原因源數(shù)據(jù)是6度帶坐標(biāo)轉(zhuǎn)換代碼里卻按3度帶中央子午線處理或者反過來。帶號與中央子午線的換算關(guān)系被忽略導(dǎo)致投影帶錯(cuò)位。解決先用已知點(diǎn)位驗(yàn)證。把源數(shù)據(jù)反算經(jīng)緯度看經(jīng)度是否落在目標(biāo)帶范圍內(nèi)再檢查Y坐標(biāo)數(shù)量級正常高斯投影Y值應(yīng)在500公里上下浮動。批量作業(yè)前務(wù)必在文件頭寫明源帶號和目標(biāo)帶號。5.2 把WGS84坐標(biāo)直接當(dāng)CGCS2000入庫現(xiàn)象數(shù)據(jù)轉(zhuǎn)完以后疊到已有CGCS2000底圖上地物整體偏移20到50厘米有的區(qū)域偏移方向還不一致旋轉(zhuǎn)縮放都救不回來。原因WGS84和CGCS2000橢球參數(shù)幾乎一樣但不同批次GPS數(shù)據(jù)的框架歷元不同導(dǎo)致系統(tǒng)性偏差。直接跳過了基準(zhǔn)轉(zhuǎn)換步驟。解決抽取項(xiàng)目區(qū)內(nèi)5個(gè)以上已知CGCS2000控制點(diǎn)用第4章的方法檢測偏移量和旋轉(zhuǎn)量。若偏差在項(xiàng)目精度允差內(nèi)可直接沿用否則必須求七參數(shù)或至少做一處平移修正。切忌不做檢核就整庫導(dǎo)入。5.3 尺度參數(shù)單位寫錯(cuò)現(xiàn)象距離越遠(yuǎn)的點(diǎn)偏差越大100公里外的點(diǎn)差了1米多。原因七參數(shù)中尺度參數(shù)m的單位是ppm對應(yīng)的無量綱數(shù)例如1ppm應(yīng)寫成0.000001。有人填成了1或者1000000尺度因子放大了百萬倍。解決轉(zhuǎn)換前先量算兩個(gè)公共點(diǎn)間的距離轉(zhuǎn)換后再次量算對比距離變化。如果距離差值與尺度參數(shù)值呈明顯相關(guān)性基本就是單位問題。正規(guī)做法是寫一個(gè)距離檢核函數(shù)轉(zhuǎn)換后自動對比多組點(diǎn)對距離。5.4 七參數(shù)和四參數(shù)模型混用現(xiàn)象同一批數(shù)據(jù)用A軟件轉(zhuǎn)換和用B軟件轉(zhuǎn)換結(jié)果差幾十厘米同行之間交換參數(shù)后各自跑出來的坐標(biāo)對不上。原因A軟件用的是七參數(shù)B軟件默認(rèn)按四參數(shù)處理或者兩家用的七參數(shù)旋轉(zhuǎn)方向約定不同。參數(shù)模型和符號約定沒有統(tǒng)一是行業(yè)合作里最常見的內(nèi)耗來源。解決統(tǒng)一轉(zhuǎn)換工具和參數(shù)文件。在項(xiàng)目中明確規(guī)定源坐標(biāo)系、目標(biāo)坐標(biāo)系、轉(zhuǎn)換模型、參數(shù)格式并在交付參數(shù)表里寫明旋轉(zhuǎn)單位弧度還是角秒、尺度單位、坐標(biāo)軸順序。多一個(gè)人用就得多一份說明這個(gè)不能省。5.5 所有公共點(diǎn)都參與擬合沒有獨(dú)立檢核點(diǎn)現(xiàn)象七參數(shù)計(jì)算時(shí)報(bào)出來的殘差只有幾毫米但到了實(shí)地放樣或與已有控制網(wǎng)比對時(shí)偏差卻達(dá)到2厘米。原因所有公共點(diǎn)都參與了參數(shù)求解轉(zhuǎn)換參數(shù)為了硬湊這些點(diǎn)而被扭曲屬于典型的過擬合。擬合殘差小不代表對其他點(diǎn)位也適用。解決至少預(yù)留三分之一的公共點(diǎn)不參與求參只用來做最終檢核。檢核點(diǎn)的中誤差才是實(shí)際可用的精度指標(biāo)。沒有條件預(yù)留時(shí)也必須至少用兩個(gè)已知點(diǎn)做旁證不能只看擬合殘差。6. 多源數(shù)據(jù)批轉(zhuǎn)與驗(yàn)證讓坐標(biāo)轉(zhuǎn)換成果經(jīng)得起復(fù)核6.1 三步驗(yàn)收法把誤差鎖死在毫米級第一步用預(yù)留檢核點(diǎn)驗(yàn)證。轉(zhuǎn)換后讀取檢核點(diǎn)的目標(biāo)坐標(biāo)與已知值逐點(diǎn)對比計(jì)算平面中誤差和最大誤差。第二步用實(shí)測RTK數(shù)據(jù)驗(yàn)證。在項(xiàng)目區(qū)選取10個(gè)以上可到達(dá)的地面點(diǎn)實(shí)測CGCS2000坐標(biāo)與轉(zhuǎn)換結(jié)果對比這一步能發(fā)現(xiàn)控制點(diǎn)坐標(biāo)本身是否存在粗差。第三步做圖面疊合驗(yàn)證把轉(zhuǎn)換成果與已有正射影像或地形圖疊加檢查道路、水系等地物是否出現(xiàn)系統(tǒng)性偏移。6.2 用NumPy向量化批量處理大CSV幾十萬條坐標(biāo)記錄逐行轉(zhuǎn)換會很慢直接把經(jīng)緯度數(shù)組傳給Transformer即可import numpy as np from pyproj import Transformer trans Transformer.from_crs( EPSG:4326, projtmerc lat_00 lon_0117 k1 x_0500000 ellpsGRS80 unitsm no_defs, always_xyTrue, ) lon np.loadtxt(lon.txt) lat np.loadtxt(lat.txt) x, y trans.transform(lon, lat) np.savetxt(xy_2000.txt, np.column_stack([x, y]), fmt%.3f)這段代碼一次轉(zhuǎn)換整批數(shù)組速度比循環(huán)快一個(gè)量級。內(nèi)存有限時(shí)可以用np.array_split分塊處理每塊10萬條轉(zhuǎn)換后寫回文件。6.3 給轉(zhuǎn)換結(jié)果留下可追溯信息坐標(biāo)轉(zhuǎn)換最大的坑不是算不出來而是算完以后說不清參數(shù)是什么。我每次交付轉(zhuǎn)換成果都會在輸出文件里附帶一列source_crs和transform_params記錄原始坐標(biāo)系、投影帶、七參數(shù)來源和版本。原始數(shù)據(jù)永遠(yuǎn)保留一份轉(zhuǎn)換結(jié)果只作為派生數(shù)據(jù)使用。這個(gè)習(xí)慣救過我很多次因?yàn)橐坏﹨?shù)更新或發(fā)現(xiàn)錯(cuò)誤還能回到原點(diǎn)重新處理。真正跑通一個(gè)轉(zhuǎn)換器靠的不是某個(gè)神秘軟件而是把橢球、投影、參數(shù)、檢核這四件事串起來。坐標(biāo)轉(zhuǎn)換這行沒有一勞永逸的萬能參數(shù)每次接手新區(qū)域都要用公共點(diǎn)重新驗(yàn)證。希望我的這些踩坑經(jīng)驗(yàn)?zāi)軒湍阍谙乱粋€(gè)數(shù)據(jù)匯交項(xiàng)目里少走彎路。本文還有配套的精品資源點(diǎn)擊獲取