用Python实现ASCII类型栅格数据转NC文件

Hello！Hello！大家好！
今天给大家带来的主要是如何利用Python读取ASCII码形式的数据，并将其转化为NC文件
预告：下一期将给大家带来ASCII码形式数据转TIF(有地理信息坐标的哦！)

数据准备

$如果无需复现程序，可以直接跳过这一部分的内容，从下一节开始啦$

为了方便本次教程的开展，我还是使用NCEP的多气压层气温数据，只是我这里为方便大家复现整个流程，现将其转为ASCII形式再用Python对其进行处理，下载网址在这里：点我，点我，快点我。

预处理

ArcGIS打开NC

接下来使用这个工具就行将TIF转为ASCII

示例数据展示

如果大家需要处理的数据是这个类型的，或者和这样子很像的文本数据（知道坐标网格的）其实处理起来的本质是一样的，那么接下里就开始本次教程把！

代码讲解

库函数准备

想要完成这个流程需要的库函数其实并不多，就下面这几个。

import xarray as xr
import numpy as np

文件读取

很简单，

def read_ASCII(file):
    """
    读取ASCII文件
    :param file: 输入是f.readlines(),一个列表文件
    :return: 还是一个列表文件   
    :功能：防止文件中间或者末尾有空行
    """

    file_tep = []
    for j in range(len(file)):
        if len(file[j].split()) != 0:
            file_tep.append(file[j])
    file = file_tep
f = open(r'D:\优快云/ascii示例数据.txt')
file = read_ASCII(f.readlines())

经纬度生成

根据ASCII文件的格式类型，我们发现前六行都是网格信息，那么我们第一步就是先要将网格信息读取出来。

for num,value in enumerate(file):
        value = value.split()
        if len(value) != 0:
            if num == 0:
                cols = int(value[-1]) #列数
            elif num == 1:
                rows = int(value[-1]) #行数
            elif num == 2:
                s_lon = float(value[-1]) #左上角经度
            elif num == 3:
                s_lat = float(value[-1]) #左上角纬度
            elif num == 4:
                cell_size = float(value[-1]) #单元格大小
            elif num == 5:
                mask_value = float(value[-1]) #掩膜数值
print(cols,rows,s_lon,s_lat,cell_size,mask_value)
lat = np.arange(-90,90+0.1,2.5)[::-1]
lon = np.arange(0,360,2.5)
print(len(lat),len(lon))
print(lat,lon)

输出结果
结果显示的是起始经纬度分别是$-1.25和-91.25$，而我在上边使用的$0和-90$，这个起始是ArcGIS本身储存网格数据的机制造成的，我们在处理的时候一般是指网格中心经纬度，而ArcGIS是网格边缘经纬度，所以就造成了相差1.25（刚好是一半的单元格大小哦！），这个不是很重要的。

数组生成

values_array = np.zeros((rows,cols)) #创建一个空矩阵
f.seek(0) # 将文件指针重新定位到文件开头
n = 0
for num,value in enumerate(file):
        value = value.split()
        if len(value) != 0: # 如果文件中有空行，则不进行读取
            if num > 5:
                for j in range(len(value)):
                    if float(value[j]) == mask_value:
                        value[j] = np.nan # 根据mask_value将mask_value设置为nan
                    else:
                        value[j] = float(value[j])
                values_array[n,:] = value
                n = n+1
    
f.close()
print(values_array)

保存为NC文件

from metpy.units import units
ds = xr.Dataset()
# *units.degC 是给气温赋单位，有没有单位都可以，锦上添花的东西
ds['1000hPa气温'] = (('lat','lon'),values_array*units.degC)
ds.coords['lat'] = ('lat',lat)
ds.coords['lon'] = ('lon',lon)
ds.to_netcdf(r'D:\优快云/ascii示例数据.nc')

完整代码奉上

import xarray as xr
import numpy as np
from metpy.units import units

def read_ASCII(file):
    """
    读取ASCII文件
    :param file: 输入是f.readlines(),一个列表文件
    :return: 还是一个列表文件   
    :功能：防止文件中间或者末尾有空行
    """

    file_tep = []
    for j in range(len(file)):
        if len(file[j].split()) != 0:
            file_tep.append(file[j])
    file = file_tep
    return file
f = open(r'D:\优快云/ascii示例数据.txt')
file = read_ASCII(f.readlines())

for num,value in enumerate(file):
        value = value.split()
        if len(value) != 0:
            if num == 0:
                cols = int(value[-1])
            elif num == 1:
                rows = int(value[-1])
            elif num == 2:
                s_lon = float(value[-1])
            elif num == 3:
                s_lat = float(value[-1])
            elif num == 4:
                cell_size = float(value[-1])
            elif num == 5:
                mask_value = float(value[-1])
print(cols,rows,s_lon,s_lat,cell_size,mask_value)
lat = np.arange(-90,90+0.1,2.5)[::-1]
lon = np.arange(0,360,2.5)
print(len(lat),len(lon))
print(lat,lon)

values_array = np.zeros((rows,cols)) #创建一个空矩阵
f.seek(0) # 将文件指针重新定位到文件开头
n = 0
for num,value in enumerate(file):
        value = value.split()
        if len(value) != 0: # 如果文件中有空行，则不进行读取
            if num > 5:
                for j in range(len(value)):
                    if float(value[j]) == mask_value:
                        value[j] = np.nan # 根据mask_value将mask_value设置为nan
                    else:
                        value[j] = float(value[j])
                values_array[n,:] = value
                n = n+1
    
f.close()
print(values_array)

ds = xr.Dataset()
# *units.degC 是给气温赋单位，有没有单位都可以，锦上添花的东西
ds['1000hPa气温'] = (('lat','lon'),values_array*units.degC)
ds.coords['lat'] = ('lat',lat)
ds.coords['lon'] = ('lon',lon)
ds.to_netcdf(r'D:\优快云/ascii示例数据.nc')

预告

下一期打算做一下如何直接使用ASCII文件生成TIF文件