C-SQUARES

C-SQUARES.md

my attempt to implement C-SQUARES in some modern languages

specs: https://www.cmar.csiro.au/csquares/spec1-1.htm

csquares.go

package main

import (
	"fmt"
	"math"
	"math/rand"
	"strconv"
	"strings"
)

func main() {
	const N = 10
	resol := [3]float64{.5, .1, .01} // degrees
	for i := 1; i <= N; i++ {
		ϕ := 180 * rand.Float64() -  90 // latitude
		λ := 360 * rand.Float64() - 180 // longitude
		fmt.Println("ϕ = ", ϕ, ", λ = ", λ)
		for _, r := range resol {
			fmt.Println(identify_csquares(ϕ, λ, r))
		}
		fmt.Println()
	}
}

/*
inputs:
  - latitude, longitude: in decimal and CRS = WGS84
  - resolution: decimal degrees, must be a number from the sequence 10, 5, 1, 0.5, 0.1, 0.05, etc.
output: C-SQUARES notation
*/
func identify_csquares(latitude, longitude, resolution float64) string {
	var (
		res strings.Builder
		i, j, k int
		ϕ_rem, λ_rem string
		ϕ_digit, λ_digit byte
	)
	k = _encode_resolution(resolution)
	// resolution=10 ⇒ k=0, resol=5 ⇒ k=1, resol=1 ⇒ k=2, resol=.5 ⇒ k=3, resol=.1 ⇒ k=4, resol=.01 ⇒ k=5, etc.
	if k == -1 {
		panic("resolution must be a number from the sequence 10, 5, 1, 0.5, 0.1, 0.05, etc.")
	}
	if k > 16 {
		panic("max precision reached with 7 decimal places in GPS coordinates (~1 cm precision)")
	}
	if latitude < -90 || latitude > 90 || longitude < -180 || longitude > 180 {
		panic("invalid GPS coordinates")
	}
	res.WriteString(_identify_quadrant(latitude, longitude, false))

	// rem = remainder
	ϕ_rem = fmt.Sprintf("%010.7f", math.Abs( latitude))
	λ_rem = fmt.Sprintf("%011.7f", math.Abs(longitude))
	// 7 decimal places in GPS coordinates for ~1 cm precision
	res.WriteString(ϕ_rem[:1] + λ_rem[:2]) // ϕ ∈ 0:8 & λ ∈ 0:17

	// make it easier for later in loop
	ϕ_rem = strings.ReplaceAll(ϕ_rem[1:], ".", "") // remove decimal separator
	λ_rem = strings.ReplaceAll(λ_rem[2:], ".", "") // make 2 strings ϕ & λ same length

	// DO NOT DO the incorrect way: separate this loop into 2 loops for even & odd index
	for i = 1; i <= k; i++ {
		j = (i + 1) / 2 - 1 // j = 0, 0, 1, 1, 2, 2, etc.
		ϕ_digit = ϕ_rem[j] // from 0 to 9
		λ_digit = λ_rem[j] // from 0 to 9
		if i % 2 == 1 { // resolution 5, 0.5, 0.05, etc.
			res.WriteString(":" + _identify_quadrant(float64(ϕ_digit), float64(λ_digit), true))
		} else { // resolution 1, 0.1, 0.01, etc.
			res.WriteString(string(ϕ_digit) + string(λ_digit))
		}
	}

	return res.String()
}

// ϕ = latitude, λ = longitude
func _identify_quadrant(ϕ, λ float64, sub_quadrant bool) string {
	var i int
	if sub_quadrant { // both ϕ & λ ∈ 0:9
		i = Btoi(λ >= 5) + 2 * Btoi(ϕ >= 5) + 1
		// 1 = south west, 2 = south east, 3 = north west, 4 = north east
	} else { // ϕ ∈ -8:8 & λ ∈ -17:17
		i = 4 * Btoi(λ < 0) + 2 * (Btoi(λ >= 0) ^ Btoi(ϕ >= 0)) + 1
		// 1 = north east, 3 = south east, 5 = south west, 7 = north west
	}
	return strconv.Itoa(i)
}

/*
let (uₖ) be the sequence u₀=10, u₁=5, u₂=1, u₃=.5, u₄=.1, u₅=.05, etc.
explicit expression of (uₖ) given by: ∀k∈ℕ
  - if k even, meaning k=2n   for some n∈ℕ, then uₖ = 10/10ⁿ
  - if k odd,  meaning k=2n+1 for some n∈ℕ, then uₖ =  5/10ⁿ
this function returns the unique integer k such that uₖ = resolution (in degree)
if no such integer exists, return -1.
*/
func _encode_resolution(degree float64) int {
	if degree <= 0 || degree > 10 {
		return -1
	}
	var x float64

	// case k even: uₖ=10/10ⁿ if k=2n
	x = math.Log10(degree) // let x=log10(uₖ) ⇒ x=1-n
	if int_check(x) {
		return 2 * (1 - int(math.Round(x))) // n=1-x ⇒ return k=2n
	}

	// case k odd: uₖ=5/10ⁿ if k=2n+1
	x -= math.Log10(5) // let x=log10(uₖ)-log10(5) ⇒ x=-n
	if int_check(x) {
		return 1 - 2 * int(math.Round(x)) // n=-x ⇒ return k=2n+1
	}

	// otherwise, no such integer k exists
	return -1
}

// check whether a float is a whole number
func int_check(x float64) bool {
	const ε = 1e-9 // margin of error
	_, x_frac := math.Modf(math.Abs(x))
	return x_frac < ε || x_frac > 1-ε
}

// golang does not have a built-in bool→int conversion
func Btoi(b bool) int {
	if b {
		return 1
	} else {
		return 0
	}
}

csquares.jl

import Printf  # do not use Printf.@sprintf coz of limited capability

""""
inputs:
- latitude, longitude: in decimal and CRS = WGS84
- resolution: decimal degrees, must be a number from the sequence 10, 5, 1, 0.5, 0.1, 0.05, etc.

output: C-SQUARES notation
"""
function identify_csquares(latitude::T, longitude::T; resolution::T)::String where {T<:Real}
	k = _encode_resolution(resolution)
	# resolution=10 ⇒ k=0, resol=5 ⇒ k=1, resol=1 ⇒ k=2, resol=.5 ⇒ k=3, resol=.1 ⇒ k=4, resol=.01 ⇒ k=5, etc.
	if k == -1 error("resolution must be a number from the sequence 10, 5, 1, 0.5, 0.1, 0.05, etc.") end
	if k > 16 error("max precision reached with 7 decimal places in GPS coordinates (~1 cm precision)") end
	if latitude ≤ -90 || latitude ≥ 90 || longitude ≤ -180 || longitude ≥ 180 error("invalid GPS coordinates") end
	res = _identify_quadrant(latitude, longitude, false)

	# rem = remainder
	ϕ_rem = Printf.format("%010.7f", abs( latitude))
	λ_rem = Printf.format("%011.7f", abs(longitude))
	# 7 decimal places in GPS coordinates for ~1 cm precision
	res *= first(ϕ_rem, 1) * first(λ_rem, 2)  # ϕ ∈ 0:8 & λ ∈ 0:17
	if k == 0 return res end

	# make it easier for later in loop
	ϕ_rem = replace(chop(ϕ_rem, head=1, tail=0), "." => "")  # remove decimal separator
	λ_rem = replace(chop(λ_rem, head=2, tail=0), "." => "")  # make 2 strings ϕ & λ same length

	# DO NOT DO the incorrect way: separate this loop into 2 loops for even & odd index
	for i ∈ 1:k
		j = (i+1) ÷ 2  # same as d above
		ϕ_digit = ϕ_rem[j]  # from 0 to 9
		λ_digit = λ_rem[j]  # from 0 to 9
		if Bool(i % 2)  # resolution 5, 0.5, 0.05, etc.
			res *= ":" * _identify_quadrant(parse(Int, ϕ_digit), parse(Int, λ_digit), true)
		else  # resolution 1, 0.1, 0.01, etc.
			res *= ϕ_digit * λ_digit
		end
	end

	return res
end


""""ϕ = latitude, λ = longitude"""
function _identify_quadrant(ϕ::T, λ::T, sub_quadrant::Bool)::String where {T<:Real}
	if sub_quadrant  # both ϕ & λ ∈ 0:9

		#= # naive version
		    if λ < 5 && ϕ ≥ 5 return "3"  # north west
		elseif λ ≥ 5 && ϕ ≥ 5 return "4"  # north east
		elseif λ ≥ 5 && ϕ < 5 return "2"  # south east
		elseif λ < 5 && ϕ < 5 return "1"  # south west
		else error("oh no")
		end
		=#

		# boolean version
		return string((λ ≥ 5) + 2*(ϕ ≥ 5) + 1)

	else  # ϕ ∈ -8:8 & λ ∈ -17:17

		#= # naive version
		    if λ < 0 && ϕ ≥ 0 return "7"  # north west
		elseif λ ≥ 0 && ϕ ≥ 0 return "1"  # north east
		elseif λ ≥ 0 && ϕ < 0 return "3"  # south east
		elseif λ < 0 && ϕ < 0 return "5"  # south west
		else error("oh no")
		end
		=#

		# boolean version
		return string(4*(λ < 0) + 2*((λ ≥ 0) ⊻ (ϕ ≥ 0)) + 1)

	end
end


@doc raw"""
let (uₖ) be the sequence u₀=10, u₁=5, u₂=1, u₃=.5, u₄=.1, u₅=.05, etc.

explicit expression of (uₖ) given by:
```math
\begin{equation}
	\forall k \in \mathbb{N},
	\begin{cases}
		u_k = \frac{10}{10^n} & \text{if k even: } \exists n \in \mathbb{N},\ k=2n\\
		u_k = \frac{ 5}{10^n} & \text{if k odd: }  \exists n \in \mathbb{N},\ k=2n+1
	\end{cases}
\end{equation}
```

this julia function returns the value of k
- input: resolution of C-SQUARES in decimal degree
- output: index k so that uₖ=input, -1 if invalid input
"""
function _encode_resolution(degree::T)::Int where {T<:Real}
	if degree < 0 || degree > 10 return -1 end

	# case k even: uₖ=10/10ⁿ if k=2n
	x = log10(degree)  # let x=log10(uₖ) ⇒ x=1-n
	try
		return 2 * (1 - Int(x))  # n=1-x ⇒ return k=2n
	catch e  # Int(x) throw InexactError if decimal part ≉ 0
		if !isa(e, InexactError) throw(e) end
	end

	# case k odd: uₖ=5/10ⁿ if k=2n+1
	x -= log10(5)  # let x=log10(uₖ)-log10(5) ⇒ x=-n
	try
		return 1 - 2 * Int(x)  # n=-x ⇒ return k=2n+1
	catch e
		if !isa(e, InexactError) throw(e) end
	end

	# otherwise, no such integer k exists
	return -1
end


if abspath(PROGRAM_FILE) == @__FILE__
	const N = 50
	const LATITUDES  = 180 .* rand(N) .- 90
	const LONGITUDES = 360 .* rand(N) .- 180
	println("ϕ = ", LATITUDES)
	println("λ = ", LONGITUDES)
	for resol ∈ [.1, .05, .01]
		println("\nresolution $resol: ", identify_csquares.(LATITUDES, LONGITUDES; resolution=resol))
	end
end