利用PYTHON将netcdf（.nc）格式数据转换成JSON GeoJSON

PYTHON netcdf（.nc）格式数据转JSON GeoJSON

• 代码
• 代码解释

代码

废话不多，直接上代码

#%%
import xarray as xr
import json
from datetime import datetime,timedelta
import numpy as np
from glob import glob

# %%
def gen_data_dict(data_array, reftime, parameterNumber, parameterNumberName):
    """
        data_array: an xarray data with two dimensions;
        parameterNumber: the number the same as example data;
        parameterNumberName: the same as example data.
            
    """
    numberPoints = np.size(data_array)
    nx, ny = data_array.shape[:]
    parameterNumberName = 'U_component_of_current'
    parameterNumber = 2
    lo1 = data_array['lon_uv'].min()
    la1 = data_array['lat_uv'].max()
    lo2 = data_array['lon_uv'].max()
    la2 = data_array['lat_uv'].min()
    dx = np.gradient(data_array['lon_uv']).mean()
    dy = np.gradient(data_array['lat_uv']).mean()
    header_dict = {
        'discipline': 10,
        'disciplineName': 'Oceanographic_products',
        'center': 0,
        'centerName': 'Ocean Modeling and Observation Laboratory',
        'refTime': reftime,
        'significanceOfRT': 0,
        'significanceOfRTName': 'Analysis',
        'parameterCategory': 1,
        'parameterCategoryName': 'Currents',
        'parameterNumber': parameterNumber,
        'parameterNumberName': parameterNumberName,
        'parameterUnit': 'm.s-1',
        'forecastTime': 0,
        'surface1Type': 160,
        'surface1TypeName': 'Depth below sea level',
        'surface1Value': 15,
        'numberPoints': numberPoints,
        'shape': 0,
        'shapeName': 'Earth spherical with radius = 6,367,470 m',
        'scanMode': 0,
        'nx': nx,
        'ny': ny,
        'lo1': float(lo1),
        'la1': float(la1),
        'lo2': float(lo2),
        'la2': float(la2),
        'dx': dx,
        'dy': dy
    }
    nan_to_none = np.fliplr(np.where(np.isnan(data_array), None, data_array))
    data_list = list(nan_to_none.ravel('F'))
    return {'header':header_dict,'data':data_list}

#%% 
if __name__ == '__main__':
    file_path = '20090101-20090103\*.nc'
    file_list = glob(file_path)
    file_list.sort()
    for file_name in file_list:
        ds = xr.open_dataset(file_name,decode_times=False)
        #%% -----------------covert time --------------
        matlab_datenum = ds['time'].values[0]
        # file_time = datetime.fromordinal(int(matlab_datenum)) + timedelta(days=matlab_datenum%1) - timedelta(days = 366)
        file_time = datetime.fromordinal(int(matlab_datenum)) - timedelta(days = 366)
        print(file_time)
        # ds['time'] = file_time
        reftime = file_time.strftime('%Y-%m-%dT00:00:00.000Z')
        # step = 10
        depth_layer = 14
        u = np.squeeze(ds['u'][depth_layer])
        parameterNumberName = 'U_component_of_current'
        parameterNumber = 2
        json_u = gen_data_dict(u, reftime,parameterNumber, parameterNumberName)
        v = np.squeeze(ds['v'][depth_layer])
        parameterNumberName = 'V_component_of_current'
        parameterNumber = 3
        json_v = gen_data_dict(v, reftime,parameterNumber, parameterNumberName)
        json_list = [json_u, json_v]
        date_str = file_time.strftime('%Y%m%d')
        if depth_layer==0:
            surface = 'surface'
        else:
            # surface = f'{depth_layer}m'
            surface = '100m'
        #%% ------------------output-------------------
        with open(f'public\data\oscar\{date_str}-{surface}-currents-oscar-0.33.json', 'w') as fp:
            json.dump(json_list, fp)
        # ds2 = xr.combine_by_coords([u,v])
        # ds2.to_netcdf(f'20090101-20090103\ROMS_0_{date_str}.nc')
    print('Done')

代码解释

首先是定义了一个函数，要求输入的第一个参数为一个二维的矩阵，实际上是个xarray的dataarray，对xarray不熟悉的可以参阅 官方文档 ，第二和第三个参数是按照给的geojson例子设定的，在本例子中u v方向分量的number分别是2和3，第三个参数是变量描述名。函数的功能是计算nx, ny, dx, dy和转换数据等功能，注意使用ravel转换成一维数据时要使用参数’F’, 才是符合JS读取的顺序的。

在实际运行中，利用xarray把数据从文件夹里挨个读取出来，然后把nc数据的时间转换一下，如果读取的时间是cftime的格式则不需要转换，可以直接读取得到datetime格式的数据，然后提取某一深度的进行转换，代码中给的是第14层，也可以加多一个对深度的循环进行循环读取，随后调用函数即可。单个JSON数据的格式其实是：

{'header':header_dict,'data':data_list}

但因为要把u v分量同时写进一个JSON中，所以则需要两个字典组成一个列表，因此使用了

json_u = gen_data_dict(u, reftime,parameterNumber, parameterNumberName)
json_v = gen_data_dict(v, reftime,parameterNumber, parameterNumberName)
json_list = [json_u, json_v]

最后利用json.dump的方式写入文件就完成了转换。