GIS 相关

PROJCS["CGCS2000 / 3-degree Gauss-Kruger CM 105E",GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",1],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","4544"]]


名称 "CGCS2000 / 3-degree Gauss-Kruger CM 105E"
表示这是以 CGCS2000 大地坐标系（中国大地2000坐标系）为基础，采用 3度带高斯-克吕格投影（中央经线105°E） 的投影坐标系（ 高斯-克吕格投影是一种将地球椭球面上的地理坐标（经纬度）正形投影到二维平面上的方法）。

GEOGCS[“China Geodetic Coordinate System 2000”, …]
表明其大地基准是 CGCS2000，是中国国家级大地坐标系，跟 WGS84 接近。


DATUM[“China_2000”, SPHEROID[“CGCS2000”,6378137,298.257222101]]
基准面采用椭球体参数：长半轴 = 6378137 m，扁率倒数 = 298.257222101。这是标准的 CGCS2000 椭球。


PROJECTION[“Transverse_Mercator”]
横轴墨卡托投影（高斯-克吕格就是横轴墨卡托的一个实现）。

latitude_of_origin = 0 → 投影起始纬度为赤道
central_meridian = 105 → 中央经线是东经 105°
scale_factor = 1 → 投影比例因子 1.0
false_easting = 500000 → 假东距 500km（防止出现负坐标）
false_northing = 0 → 假北距 0

UNIT[“metre”,1]
投影结果的单位是米。

AUTHORITY[“EPSG”,“4544”]
EPSG 编号 4544，对应 CGCS2000 / 3-degree Gauss-Kruger zone 35（中央经线105°E）


以上信息说明：EPSG:4544 投影坐标系的“原 点”对应 (105°E, 0°N)，中央子午线与赤道的交点，X轴：中央子午线的投影。Y轴：赤道的投影。
但因为约定俗称，要有 false easting = 500,000 m，所以这个点在平面上是 (500000, 0)。

在这个投影坐标系上，正射影像上每个像素的坐标确实是以 米为单位的 (E,N)。所以，如果输出的位置给了 (E、N)这样的坐标，需要转换成 EPSG:4490 才能用于地图显示。使用 gdal 工具：

# 使用 gdaltransform（把 E N 放在 stdin）
echo "605732.799 3999843.061 1623.661" | gdaltransform -s_srs EPSG:4544 -t_srs EPSG:4490
# 通常输出：<lon> <lat> <h>

$ gdalinfo  dzh.tif

Driver: GTiff/GeoTIFF
Files: dzh.tif
Size is 86842, 67258
Coordinate System is:
PROJCRS["CGCS2000 / 3-degree Gauss-Kruger CM 105E",
    BASEGEOGCRS["China Geodetic Coordinate System 2000",
        DATUM["China 2000",
            ELLIPSOID["CGCS2000",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4490]],
    CONVERSION["Gauss-Kruger CM 105E",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",105,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],
        AREA["China - between 103°30'E and 106°30'E."],
        BBOX[22.5,103.5,42.21,106.51]],
    ID["EPSG",4544]]
Data axis to CRS axis mapping: 2,1
Origin = (602832.567790000000969,4001299.158060000278056)
Pixel Size = (0.037220000000000,-0.037220000000000)
Metadata:
  TIFFTAG_SOFTWARE=pix4dmapper
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
  PREDICTOR=2
Corner Coordinates:
Upper Left  (  602832.568, 4001299.158) (106d 8'32.85"E, 36d 8'11.59"N)
Lower Left  (  602832.568, 3998795.815) (106d 8'31.68"E, 36d 6'50.39"N)
Upper Right (  606064.827, 4001299.158) (106d10'42.11"E, 36d 8'10.34"N)
Lower Right (  606064.827, 3998795.815) (106d10'40.89"E, 36d 6'49.13"N)
Center      (  604448.697, 4000047.487) (106d 9'36.88"E, 36d 7'30.37"N)
Band 1 Block=86842x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=86842x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=86842x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=86842x1 Type=Byte, ColorInterp=Alpha

gdalwarp -t_srs EPSG:3857 \
  -tr 1 1 \
  -r bilinear \
  -of GTiff \
  -co COMPRESS=DEFLATE \
  -co TILED=YES \
  -dstalpha \
  input.tif output_3857.tif

gdal2tiles.py -z 0-18 -r bilinear -p mercator -a 0 \
  --processes=4 warped.tif tiles/ --tilesize=256 --tiledriver=PNG

import math
from pyproj import Proj, Transformer

# 输入：两端点经纬度（度）
# lat1, lon1 = 36.134616757252154,106.13848686218262
# lat2, lon2 = 36.13080407775483,106.18852615356447

lat1, lon1 = 36.13822130205884,106.17324829101564
lat2, lon2 = 36.11617553449072,106.17470741271974

# 定义 CGCS2000 / 3-degree Gauss-Kruger 投影
# 这里假设中央子午线为 105E
proj = Proj(proj='tmerc', lat_0=0, lon_0=105, k=1, x_0=500000, y_0=0, ellps='GRS80')

# 经纬度转平面坐标
x1, y1 = proj(lon1, lat1)
x2, y2 = proj(lon2, lat2)

# 坝轴向量
dx = x2 - x1
dy = y2 - y1

# 计算方位角 θ（顺时针从正北算起）
theta_rad = math.atan2(dx, dy)  # 注意 dx 对应东向，dy 对应北向
theta_deg = math.degrees(theta_rad)
if theta_deg < 0:
    theta_deg += 360

print(f"坝轴方位角 θ = {theta_deg:.3f}°")%

// GNSS设备上报的数据。
// DzhGnssRealtimeData 保存dzhGnss设备的最新实时值
type DzhGnssRealtimeData struct {
	Lat       float64 `json:"lat"`       // 纬度
	Lon       float64 `json:"lon"`       // 经度
	Alt       float64 `json:"alt"`       // 海拔高度，单位 m
	X         float64 `json:"x"`         // X方向位移，单位 mm
	Y         float64 `json:"y"`         // Y方向位移，单位 mm
	Z         float64 `json:"z"`         // Z方向位移，单位 mm
	Power     float64 `json:"power"`     // 电压
	RSSI      int     `json:"rssi"`      // 信号强度
	Timestamp int64   `json:"timestamp"` // Unix时间戳
}

/*
"data":"{
	"lat":36.123160717,
	"lon":106.17444795,
	"alt":1620.93103,
	"drft":0.0047,
	"drfH":270,
	"drfV":0,
	"mod":4,
	"pow":12.3,
	"RSSI":22,
	"clientId":"sv78dRYjrD1eICq51QcD",
	"detectedTime":"2025-08-11 13:31:32"
}
drft 是位移大小，单位是米（m），两个点之间的间距
drfH 是水平角，单位是度（°），这里假设角度是以正北为0度，顺时针方向
drfV 是垂直角，单位是度（°），正值表示向上仰角，负值表示向下

drft = 位移大小（径向距离 r）
drfH = 水平方位角（θ，北 0° 顺时针）
drfV = 垂直仰角（φ，水平 0°，上正下负）

位移量，单位米
水平角度，单位度，此角度为设备所在水平面的位移方向，正北为0度，按顺时针增大，正东为90度
垂直角度，单位度，此角度为设备所在垂直方向的位移方向，水平为0度，向下最大垂直方向-90度，向上最大垂直方向90度
*/


// 位移(大小+水平角+垂直角) 转换为 ENU 三维向量
// drft 单位 m，drfH 水平角(°，北0°顺时针)，drfV 垂直角(°，正上仰角)
// （北+、东+、上+），相反是负的

// 球面位移转 ENU
func PolarToVector(drft, drfH, drfV float64) (dN, dE, dU float64) {
	// 角度转弧度
	hRad := drfH * math.Pi / 180.0
	vRad := drfV * math.Pi / 180.0

	// 分解
	hLen := drft * math.Cos(vRad) // 水平投影
	dE = hLen * math.Sin(hRad)    // 东向   	E (East)：向东（经度增大的方向）
	dN = hLen * math.Cos(hRad)    // 北向       N (North)：向北（纬度增大的方向）
	dU = drft * math.Sin(vRad)    // 垂直分量    U (Up)：向上（椭球面法线方向，即大地高程方向）

	// x, y z := dE,dN, dU

	return
}


// ENU 转大坝坐标系
// Y轴：沿坝轴线方向，面朝上游，右侧为正，左侧为负
// X轴：垂直于坝轴线，下游为正，上游为负
// Z轴：竖直向上

func VectorToDam(dN, dE, dU, theta float64) (dX, dY, dZ float64) {
	// θ = 坝轴线相对正北的方位角，顺时针为正
	thRad := theta * math.Pi / 180.0

	// 沿坝轴方向 → Y
	dY = dN*math.Cos(thRad) + dE*math.Sin(thRad)

	// 垂直坝轴方向 → X（下游为正）
	dX = -dN*math.Sin(thRad) + dE*math.Cos(thRad)

	// 垂直方向保持不变
	dZ = dU

	return
}