prepare for adw

src/Interpolation/angle.jl

@@ -0,0 +1,40 @@
+export angle_azimuth, angle_azimuth_sphere
+
+
+"""
+球面方位角 (0~360°)
+
+```julia
+angle_azimuth_sphere((0, 0), (0, 1)) # 正北，0
+angle_azimuth((0, 0), (0, 1))        # 正北，0
+```
+"""
+function angle_azimuth_sphere(p1::Tuple{FT,FT}, p2::Tuple{FT,FT}) where {FT}
+  lon1, lat1 = p1
+  lon2, lat2 = p2
+
+  φ1 = deg2rad(lat1)
+  φ2 = deg2rad(lat2)
+  Δλ = deg2rad(lon2 - lon1)
+
+  y = sin(Δλ) * cos(φ2)
+  x = cos(φ1) * sin(φ2) - sin(φ1) * cos(φ2) * cos(Δλ)
+  θ = atan(y, x)
+
+  azimuth = rad2deg(θ)
+  azimuth < 0 && (azimuth += 360) # 标准化方位角为0-360度
+  return azimuth
+end
+
+function angle_azimuth(p1::Tuple{FT,FT}, p2::Tuple{FT,FT}) where {FT}
+  lon1, lat1 = p1
+  lon2, lat2 = p2
+  dx = lon2 - lon1
+  dy = lat2 - lat1
+
+  # 需要atan2(Δx, Δy)来计算从北向东的方位角
+  # azimuth = -rad2deg(atan(dy, dx)) + 90
+  azimuth = rad2deg(atan(dx, dy))
+  azimuth < 0 && (azimuth += 360) # 标准化方位角为0-360度
+  return azimuth
+end

src/Interpolation/utilize.jl

@@ -0,0 +1,23 @@
+"""
+rdist_earth(x1::Matrix, x2::Matrix; R=6378.388)
+
+```julia
+p1 = [110 30; 111 31]
+p2 = [113 32; 115 35]
+rdist_earth(p1, p2)
+```
+"""
+function rdist_earth(x1::Matrix, x2::Matrix; R=6378.388)
+  lon1 = deg2rad.(x1[:, 1])
+  lat1 = deg2rad.(x1[:, 2])
+  lon2 = deg2rad.(x2[:, 1])
+  lat2 = deg2rad.(x2[:, 2])
+
+  _x = @. [cos(lat1) * cos(lon1) cos(lat1) * sin(lon1) sin(lat1)]
+  _y = @. [cos(lat2) * cos(lon2) cos(lat2) * sin(lon2) sin(lat2)]
+  pp = _x * _y'
+  return R .* acos.(clamp.(pp, -1, 1))
+end
+
+
+export rdist_earth

test/interp/test-angle.jl

@@ -0,0 +1,13 @@
+using Test
+
+@testset "azimuth" begin
+  @test angle_azimuth_sphere((0, 0), (0, 1)) == 0    # 正北，0
+  @test angle_azimuth_sphere((0, 0), (1, 0)) == 90   # 正东，90
+  @test angle_azimuth_sphere((0, 0), (0, -1)) == 180 # 正南，180
+  @test angle_azimuth_sphere((0, 0), (-1, 0)) == 270 # 正西，270
+
+  @test angle_azimuth((0, 0), (0, 1)) == 0    # 正北，0
+  @test angle_azimuth((0, 0), (1, 0)) == 90   # 正东，90
+  @test angle_azimuth((0, 0), (0, -1)) == 180 # 正南，180
+  @test angle_azimuth((0, 0), (-1, 0)) == 270 # 正西，270
+end

test/runtests.jl

@@ -6,12 +6,13 @@ dir_root = dirname(dirname(@__FILE__))
 proj_path(f) = dirname(dirname(@__FILE__)) * "/" *  f
 # println(dirname(@__FILE__))
 # println(pwd())
+include("interp/test-bilinear.jl")
+include("interp/test-angle.jl")
 
 # cd(dirname(@__FILE__)) do
 include("test-nc_date.jl")
 include("test-utilize.jl")
 include("test-ncatt.jl")
-include("test-bilinear.jl")
 include("test-MFDataset.jl")
 include("test-nc_write.jl")
 include("test-nc_dims.jl")