# SRTMReader

**Repository Path**: JetLab/srtmreader

## Basic Information

- **Project Name**: SRTMReader
- **Description**: 使用C#读取SRTM高程数据文件hgt, 并生成灰度图
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2023-12-11
- **Last Updated**: 2023-12-11

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# NASA SRTM 高程数据C#代码读取实现

<!-- TOC -->

- [NASA SRTM 高程数据C#代码读取实现](#nasa-srtm-高程数据c代码读取实现)
  - [1. 概述](#1-概述)
  - [2. SRTM](#2-srtm)
    - [2.1. 基本特性](#21-基本特性)
    - [2.2. 文件特性](#22-文件特性)
    - [2.3. 数据下载](#23-数据下载)
  - [3. 直接硬读](#3-直接硬读)
  - [4. CSharp代码实现](#4-csharp代码实现)
  - [5. 测试代码](#5-测试代码)

<!-- /TOC -->

## 1. 概述

本文写作于2023年12月09日.


因最近打算学习使用高程灰度图+卫星影像照片, 直接通过godot引擎渲染器生成三维影像, 特此记录高程数据读取至生成灰度图的过程, 暂时只有8848/255≈34.7m的高度精度, 但是考虑到SRTM本身每个像素点的水平误差也为30m, 对于学习来说精度可用.

NASA官网可以同时下载到原始数据(.hgt)和灰度图文件(.jpg), 但经分析发现, 官方提供的单张图片应该是均基于本图片经纬度原始数据范围内的最小和最大高程进行了归一化, 会导致多张图片拼接时出现偏差, 直接渲染成三维必然出现裂缝.

因此, 本人选择了自己通过编写C#代码读取hgt文件生成灰度图的思路, NASA灰度图和本人生成的灰度图效果见如下2图.


![NASA原图](image/readme/NASA原图.png)


![自编代码拼图效果](image/readme/自编代码拼图效果.png)

## 2. SRTM

[官方介绍网站](https://lpdaac.usgs.gov/products/srtmgl1v003/):https://lpdaac.usgs.gov/products/srtmgl1v003/

![官方介绍](image/readme/官方介绍.png)

### 2.1. 基本特性
- SRTMGL1v003(版本?)
- 每像素点1角秒(约30米)
- 每个文件覆盖范围为1°×1°
- 数据库覆盖范围北纬60°到南纬56°, 西经180°到东经180°

### 2.2. 文件特性
- 后缀为.hgt
- 文件名形式为NxxEyyy.hgt(yyy:经度以3位数表示, 不足3位补零)
- 文件名表示了高程数据的左下角(北纬， 东经上看), 右上角为文件名经纬度各＋1
- 以地球东北部分来看, 高度数据先经度增大, 后纬度增大编排, 因此打开这个二进制文件, 读取到的第1个数据点为文件名纬度+1位置的高程(后面要考)
- 数据总数为3601×3601
- 每个数据点由2个byte组成的有符号数, 且为大端编码
- 右键可以查看到文件的大小为25,934,402字节(即3601×3601×2)

- 数据有效范围-32767~32767(但是最高点珠穆朗玛峰只到8848?)

![文件属性](image/readme/文件属性.png)

### 2.3. 数据下载
[官方下载地址](https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/):https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/

但是现在网络上有按区域下载好的压缩包

## 3. 直接硬读

在有高程信息的地图工具上找一个参考点, 比如N31E104, 显示的海拔为571m  
![N31E104](image/readme/N31E104.jpg)

下载该区域的文件N30E104.hgt(为什么不是N31E104在2.2文件特性已说), 读取第1,2个字节分别为0x02 0x3C, 即572m与从地图上采到的差不多

![前2个字节](image/readme/前2个字节.png)

## 4. CSharp代码实现
```cs
using System;
using System.Diagnostics;
using System.Drawing; // 注意项目右键 > 添加> 引用, 添加这个程序集, 仅using不行
using System.IO;

namespace SRTMReader
{
    /*    
     * < hgt文件编排规则 >
     * - 一般的地图, 以地球东北部分为例, 单个瓦片向右经度变大, 向上纬度变大
     * - hgt文件是一个2字节连续编排的纯二进制文件, 第一个读到的数并不是文件名所在经纬度
     *   左上角--经度变大->右上角
     *         |            |
     *    纬度变大           |
     *         ↓            ↓
     *   左下角(文件名)----->右下角
     *
     * < 单个点读取规则 >
     * - srtm数据, 1个文件跨度为1°（经/纬都是）, 文件内包含3601×3601个点
     * - 每个点由2个字节的有符号整数编码而成
     * - 二进制数为大端编码, 即高8位在低地址
     * - 右键查看hgt文件时, 大小应该为25,934,402字节(即3601×3601×2)
     *
     * < SRTM高程图原始灰度图规则 >
     * - 从NASA网站下载到的图片为越黑, 海拔越低
     */
    public class Read
    {
        /// <summary>
        /// SRTMGL1V003的.hgt文件中行, 列数和行数都为3601
        /// </summary>
        public static int Size { get; set; } = 3601;

        /// <summary>
        /// 渲染参考的最大高度
        /// </summary>
        public static int HeightMax { get; set; } = 8925;

        /// <summary>
        /// 以二进制形式读取hgt文件到有符号16进制数据
        /// 注意原始数据是大端编码, 先读到的是高8位然后是低8位
        /// </summary>
        /// <param name="filePath">hgt文件路径</param>
        /// <param name="isNormalized">是否将当前文件数据规范化, 以本文件进行规范化, 可能高度差在灰度上对比更强烈</param>
        /// <param name="scale">比例尺, 只有指定进行规范化才有意义, 返回的是整型, 但规范化后范围理论上只能是-1~1, 用这个scale乘以后将数放大以使得整型输出可以保证足够精度</param>
        public static short[] ToInt16Array(string filePath, bool isNormalized = false, int scale = 8925)
        {
            var bytes = File.ReadAllBytes(filePath);
            
            var data = new short[bytes.Length / 2];

            if (bytes.Length != Size * Size * 2) Debugger.Break();

            if (isNormalized) // 根据本地的最大最小值进行数据规范化
            {
                double lowest = short.MaxValue;
                double highest = -32767;
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = (short)(bytes[i * 2] << 8 | bytes[i * 2 + 1]);
                    if (lowest > data[i])
                        lowest = data[i];

                    if (highest < data[i])
                        highest = data[i];
                }

                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = (short)(((data[i] - lowest) / (highest - lowest)) * scale);
                }
            }
            else
            {
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = (short)(bytes[i * 2] << 8 | bytes[i * 2 + 1]);
                }
            }

            return data;
        }

        /// <summary>
        /// 保存为灰度图的实现, 目前数据范围处理到0~8925m
        /// 以8位深度颜色保存则比例尺为8925/255=35m
        /// TODO: 优化范围
        /// TODO: 实现负海拔渲染
        /// TODO: 以16位以上位宽显示的实现
        /// </summary>
        /// <param name="filePath">hgt文件路径</param>
        /// <param name="outFilePath">输出的灰度图路径</param>
        /// <param name="isNormalized">是否将当前文件数据规范化, 以输出类似NASA原始图风格的图片</param>
        public static void ToGrayScale(string filePath, string outFilePath, bool isNormalized = false)
        {
            var data = ToInt16Array(filePath, isNormalized, HeightMax);

            Color[] grayScale = new Color[data.Length];

            double val;
            int r, g, b;
            for (int i = 0; i < data.Length; i++)
            {
                val = data[i];
                if (val > 0)
                    val = (int)(Math.Min(val, HeightMax) / (HeightMax) * 255);
                else
                    val = 0; // 暂时不处理负海拔

                r = (int)(val);
                b = (int)(val);
                g = (int)(val);

                grayScale[i] = Color.FromArgb(r, g, b);
            }

            Bitmap img = new Bitmap(Size, Size);
            for (int y = 0; y < Size; y++)
            {
                for (int x = 0; x < Size; x++)
                {
                    img.SetPixel(x, y, grayScale[y * Size + x]);
                }
            }

            img.Save(outFilePath);
        }
    }
}


```

## 5. 测试代码
```cs
using Microsoft.VisualStudio.TestTools.UnitTesting;

namespace SRTMReader.Test
{
    [TestClass]
    public class ReadTest
    {
        readonly string region = "N31E103";

        [TestMethod]
        public void CanConvertToInt16()
        {
            var data = SRTMReader.Read.ToInt16Array($"{region}.hgt");

            Assert.IsNotNull(data);
            Assert.AreEqual(3601 * 3601, data.Length);
        }

        [TestMethod]
        public void CanCreateGrayScale()
        {
            string grayScaleImage = $"{region}.png";
            SRTMReader.Read.ToGrayScale($"{region}.hgt", grayScaleImage);
        }

        [TestMethod]
        public void CanCreateGrayScaleNasaLike()
        {
            string grayScaleImageNasaLike = $"{region}_nasa_like.png";
            SRTMReader.Read.ToGrayScale($"{region}.hgt", grayScaleImageNasaLike, true);
        }
    }
}


```

## 6. 生成效果

### 6.1. 以8925为高度最大值
![N31E103](image/readme/N31E103.png)
### 6.2. 以本数据文件最大、最小值为极值
![N31E103_nasa_like](image/readme/N31E103_nasa_like.png)