使用高程数据

近期有个需求，在地图上绘制高度数据，于是各种搜索，期间想到的方法主要包括以下几种：

1. 使用WRF模型的WPS直接输出高度信息，自行设定好网格和投影即可；
2. 使用下载的高程数据STRM3，可以在这里下载欧亚大陆的高程数据并进行处理获取某个经纬度的高度数据；关于STRM3的介绍可以查看度娘，有很多中文的介绍；

由于当前使用的计算机没有安装WRF且需求比较着急，因此第一条方案暂且略过；
采用第二种方案，那么打开欧亚大陆的具体数据文件地址进行下载，大约有5876个ZIP文件需下载，手动实在麻烦，而且在Win环境下不是那么方便，因此网上直接搜索wget for windows，下载所需要的wget版本并加入环境变量，或者直接将wget.exe放入系统盘的windows/system32下也是可以的；

wget -c -r -nc https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/

在终端输入这条命令即开始下载，当然这么下载是会递归所有目录的，没关系，将来直接使用；

下载下来的是zip文件还无法直接使用，于是等下载好之后批量解压吧；

解压之后使用如下网络搜索的代码:

import os
import json
import numpy as np

SAMPLES = 1201  # Change this to 3601 for SRTM1
HGTDIR = 'hgt'  # All 'hgt' files will be kept here uncompressed


def get_elevation(lon, lat):
    hgt_file = get_file_name(lon, lat)
    if hgt_file:
        return read_elevation_from_file(hgt_file, lon, lat)
    # Treat it as data void as in SRTM documentation
    # if file is absent
    return -32768


def read_elevation_from_file(hgt_file, lon, lat):
    with open(hgt_file, 'rb') as hgt_data:
        # HGT is 16bit signed integer(i2) - big endian(>)
        elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES*SAMPLES)\
                                .reshape((SAMPLES, SAMPLES))

        lat_row = int(round((lat - int(lat)) * (SAMPLES - 1), 0))
        lon_row = int(round((lon - int(lon)) * (SAMPLES - 1), 0))

        return elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)

def get_file_name(lon, lat):
    """
    Returns filename such as N27E086.hgt, concatenated
    with HGTDIR where these 'hgt' files are kept
    """

    if lat >= 0:
        ns = 'N'
    elif lat < 0:
        ns = 'S'

    if lon >= 0:
        ew = 'E'
    elif lon < 0:
        ew = 'W'

    hgt_file = "%(ns)s%(lat)02d%(ew)s%(lon)03d.hgt" % {'lat': abs(lat), 'lon': abs(lon), 'ns': ns, 'ew': ew}
    hgt_file_path = os.path.join(HGTDIR, hgt_file)
    if os.path.isfile(hgt_file_path):
        return hgt_file_path
    else:
        return None

# Mt. Everest
print get_elevation(86.925278, 27.988056)
# Kanchanjunga
print get_elevation(88.146667, 27.7025)

经过上面的一系列处理，基本上我们已经得到的大部分的信息，那么下一步就开始得到网格化的高度数据.

import numpy as np

lons = np.arange(70,140 + 1e-10, 0.01)
lats = np.arange(20, 60 + 1e-10, 0.01)
gridlon, gridlat = np.meshgrid(lons, lats)
points = np.columns_stack((gridlon.ravel(), gridlat.ravel()))
hgts = np.array([
    get_elevation(*point) for point in points
]).reshape(gridlon.shape)

import xarray as xr

da_hgt = xr.DataArray(hgts, coords=(lons, lats), dims=['lon', 'lat'])

da_hgt.to_netcdf('hgt_STRM3.ng') 

到这里我们就把结果进行了保存，由于目前还没有测试上述代码，后面继续看看与其他结果进行相应的对比。

待后续进展再进行补充。