WRS2数据查询代码(Golang)

WRS2(全球参考系统)简介

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())
		}
	}
}