WRS2數據查詢代碼(Golang)

WRS2(全球參考系統)簡介

參考：全球參考系統（WRS）概述及常用工具匯總

WRS（Worldwide Reference System）是Landsat衛星采用的全球參考系統，也是國際上非常具有代表意義的全球參考系統之一。WRS是依據衛星地面軌跡的重復特性，結合星下點成像特性而形成的固定地面參考網格。其WRS參考系網格與Landsat衛星數據的成像區域緊密的契合，WRS網格的二維坐標采用Path和Row進行標識。
目前WRS有兩個系統，分別為WRS-1（1983年之前）和WRS-2（1983年之后）。Landsat1-3使用WRS-1，Landsat4、5、7、8使用WRS-2。

WRS2 數據

這個數據我找了很久，NASA的網站現在只有在線查詢的地方，沒有找到數據下載。后來在這里 Converting latitude/longitude co-ordinates to Landsat WRS-2 paths/rows 找到了數據的下載位置。

代碼

完整代碼請見：https://github.com/sotex/Landsat8cvsToMySQL

WRS2 查詢相關代碼如下：

package main

import (
	"bytes"
	"encoding/json"
	"fmt"
	"strconv"
	"strings"

	"github.com/tidwall/geojson"

	"github.com/jonas-p/go-shp"
	geometry "github.com/tidwall/geojson/geometry"
)

type (
	WRS2Object struct {
		Path     int            // 路徑
		Row      int            // 行
		boundary *geometry.Poly // 有效范圍
	}
)

// 轉換為 geojson 格式
func (obj *WRS2Object) MarshalJSON() ([]byte, error) {
	buf := bytes.NewBuffer(nil)
	buf.WriteString(`{"type":"Feature","properties":{"path":`)
	buf.WriteString(strconv.Itoa(obj.Path))
	buf.WriteString(`,"row":`)
	buf.WriteString(strconv.Itoa(obj.Row))
	buf.WriteString(`},"geometry":`)
	data, err := json.Marshal(geojson.NewPolygon(obj.boundary))
	if err != nil {
		return nil, err
	}
	buf.Write(data)
	buf.WriteByte('}')
	return buf.Bytes(), err
}

var (
	// WRS2索引 key=path*1000+row
	sPathRowIndex map[int]*WRS2Object
)

func init() {
	sPathRowIndex = make(map[int]*WRS2Object)
}

// 加載 wrs2 數據
// 這里加載的數據是經過處理的，把極地附近的給去掉了
func loadWRS2Data(shpfile string) error {
	shape, err := shp.Open(shpfile)
	if err != nil {
		return fmt.Errorf("打開文件錯誤:%s", err.Error())
	}
	defer shape.Close()

	// 字段信息
	// AREA PERIMETER PR_ PR_ID RINGS_OK RINGS_NOK PATH ROW MODE SEQUENCE WRSPR PR ACQDayL7 ACQDayL8
	// fields := shape.Fields()
	// fmt.Println(fields)
	for shape.Next() {
		n, s := shape.Shape()
		p, ok := s.(*shp.Polygon)
		if !ok {
			continue
		}
		// 獲取邊界點坐標
		coordinates := make([][]geometry.Point, p.NumParts)
		for i := int32(0); i < p.NumParts; i += 1 {
			var startIndex, endIndex int32
			startIndex = p.Parts[i]
			if i == p.NumParts-1 {
				endIndex = int32(len(p.Points))
			} else {
				endIndex = p.Parts[i+1]
			}

			coordinates[i] = make([]geometry.Point, 0, endIndex-startIndex)
			cnt := 0
			for j := startIndex; j < endIndex; j = j + 1 {
				var point = geometry.Point{
					X: p.Points[j].X,
					Y: p.Points[j].Y,
				}
				// 如果不是第一個或者最后一個點，判斷與前一個點的距離
				// 用於減少點數（變成四邊形，雖然不准確）
				if cnt > 0 && cnt < int(endIndex-startIndex-1) {
					dx, dy := point.X-coordinates[i][cnt-1].X, point.Y-coordinates[i][cnt-1].Y
					// 如果距離小於0.1度，就認為是一個點
					if (dx*dx + dy*dy) < 0.1 {
						continue
					}
				}
				coordinates[i] = append(coordinates[i], point)
				cnt = cnt + 1
			}
		}
		// 獲取路徑和行號
		var path, row int
		path, _ = strconv.Atoi(shape.ReadAttribute(n, 6))
		row, _ = strconv.Atoi(shape.ReadAttribute(n, 7))
		// 插入map里面
		var pr = path*1000 + row
		polygon := geometry.NewPoly(coordinates[0], coordinates[1:], nil)
		sPathRowIndex[pr] = &WRS2Object{
			Path:     path,
			Row:      row,
			boundary: polygon,
		}
	}
	return nil
}

func queryWRS2(path, row int) *WRS2Object {
	p, ok := sPathRowIndex[path*1000+row]
	if ok {
		return p
	}
	return nil
}

// 簡單的測試代碼
func main() {
	loadWSR2Data("E:/WRS2_descending_0/WRS2_descending.shp")
	fmt.Println("數據加載完成,共讀取:", len(sPathRowIndex))

	p := queryWSR2(13, 10)
	if p != nil {
		data, err := json.Marshal(p)
		if err == nil {
			fmt.Println(string(data))
		} else {
			fmt.Println(err.Error())
		}
	}
}