python Xarray處理設置二維數組作為coordinates方式

更新時間：2023年07月07日 11:05:15 作者：永康李

這篇文章主要介紹了python Xarray處理設置二維數組作為coordinates方式，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教

python Xarray處理設置二維數組作為coordinates

因為想做筆記，所以直接做的很粗糙了，后面再更新！

import cv2
import numpy as np
from osgeo import gdal
import os
import xarray as xr 
import matplotlib.pyplot as plt
import matplotlib as mpl
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
cmap = mpl.cm.cool
norm = mpl.colors.Normalize(vmin=5, vmax=10)
fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),
             cax=ax, orientation='horizontal', label='Some Units')
"""
res = cv2.resize(RasterArrray, dsize=(441,251), interpolation=cv2.INTER_CUBIC)
Here img is thus a numpy array containing the original image, whereas res is a numpy array containing the resized image. An important aspect is the interpolation parameter: there are several ways how to resize an image. Especially since you scale down the image, and the size of the original image is not a multiple of the size of the resized image. Possible interpolation schemas are:
INTER_NEAREST - a nearest-neighbor interpolation
INTER_LINEAR - a bilinear interpolation (used by default)
INTER_AREA - resampling using pixel area relation. It may be a preferred method for image decimation, as it gives moire'-free results. But when the image is zoomed, it is similar to the INTER_NEAREST method.
INTER_CUBIC - a bicubic interpolation over 4x4 pixel neighborhood
INTER_LANCZOS4 - a Lanczos interpolation over 8x8 pixel neighborhood
"""
def GetTimeSerises_nc(ncVariable):
    """
    獲取 時間序列
    :param ncVariable:
    :return:
    """
timeSerises = ncVariable.time.data
return timeSerises
inNcFile = r"./solar-1979-01.nc"
inNc = xr.open_dataset(inNcFile)
print(inNc)
print(inNc.LATIXY.data)
import pandas as pd 
# 創(chuàng)建 dataset
ds = xr.Dataset()
numLon = 1400
numLat = 800
# LATIXY LONGXY
inLat = inNc.LATIXY.data
inLon = inNc.LONGXY.data
# print("np.min(inLon):{}, np.max(inLon):{}".format(np.min(inLon), np.max(inLon)))
# print("np.min(inLat):{}, np.max(inLat):{}".format(np.min(inLat), np.max(inLat)))
lon = np.linspace(np.min(inLon), np.max(inLon), num=numLon, endpoint=True, retstep=False, dtype=None, axis=0)
lat = np.linspace(np.min(inLat), np.max(inLat), num=numLat, endpoint=True, retstep=False, dtype=None, axis=0)
lon, lat = np.meshgrid(lon, lat)
ds = ds.assign_coords({
    "lat": (["x", "y"], lat),
    "lon": (["x", "y"], lon)
})
solor = np.full(shape=(10, numLat, numLon) , fill_value= np.nan )
ncVariable = inNc.FSDS
timeSerises = GetTimeSerises_nc(ncVariable)
i = 0
for timeSerise in timeSerises[0:10]:
    print(timeSerise)
    # 獲取數據
    arr = inNc.FSDS.loc[timeSerise].data
    print(arr.shape)
    solor[i,:,:] = cv2.resize(arr, dsize=(numLon,numLat), interpolation = cv2.INTER_LINEAR)
    print(arr.shape)
    i= i+1
    print(i)
ds["solor"] = xr.DataArray(solor, dims=['time','x', 'y'], )
ds.coords['time'] = pd.date_range(start='1979-01-01',periods=10,freq='3H')
# ds["lat"]  = xr.DataArray(lat, dims=['lat'], )
# ds["lon"]  = xr.DataArray(lon, dims=['lon'], )
print(ds)
ds.to_netcdf(r"./test_1.nc")

主要解決問題的代碼塊在這里：

lon = np.linspace(np.min(inLon), np.max(inLon), num=numLon, endpoint=True, retstep=False, dtype=None, axis=0)
lat = np.linspace(np.min(inLat), np.max(inLat), num=numLat, endpoint=True, retstep=False, dtype=None, axis=0)
lon, lat = np.meshgrid(lon, lat)
ds = ds.assign_coords({
    "lat": (["x", "y"], lat),
    "lon": (["x", "y"], lon)
})
ds["solor"] = xr.DataArray(solor, dims=['time','x', 'y'], )
ds.coords['time'] = pd.date_range(start='1979-01-01',periods=10,freq='3H')

結果：

參考鏈接https://stackoverflow.com/questions/67695672/xarray-set-new-2d-coordinate-as-dimension

Xarray(python)讀取?Sentinel-5P(S5P)哨兵數據

需求分析：NC文件的常規(guī)包netcdf4使用手感較xarray略顯笨拙，故嘗試使用xarray讀取包含Group的.nc4文件

數據：S5P二級數據：S5P_RPRO_L2__HCHO，來源：歐洲哥白尼，或NASA（推薦，因為好下載）

使用panoly可視化

（1）導入后的界面：

（2）選擇變量后，點擊Create Plot按鈕可視化：

即可得到HCHO的Plot圖以及Array可視化。

使用python里的工具包讀取

import os
import xarray as xr
import netCDF4 as nc  # 對于nc4文件，其內含groups，
Dir = ['../S5P_Pre/Wget_HCHO']   # 時間跨度180514 ~ 190805
file = os.listdir(Dir[0])
file.sort(key = lambda x:int(x.split('___')[1][:8]))  # 按年月日排序
# （1）使用nc包打開
ns = nc.Dataset(os.path.join(Dir[0], file[0]))   #這里的數據存儲在groups里面的PRODUCT里面
hcho = ns['PRODUCT']['formaldehyde_tropospheric_vertical_column'][:]
# (2) 使用xarray包打開 —— 推薦方式
xs = xr.open_dataset(os.path.join(Dir[0], file[0]), group = 'PRODUCT')  # 這里需用group函數指定組名稱

（1）netcdf4的讀取結果：

In[29]: ns
Out[29]: Subset parameters: {"PRODUCT": ["S5P_L2__HCHO__.1"], "INFILENAMES": ["S5P_RPRO_L2__HCHO___20180514T023918_20180514T042246_03018_01_010105_20190203T205044.nc"], "INFILETYPE": ["nc"], "OUTFILETYPE": ["nc4"], "TIMENAME": [["TROP2010", "/PRODUCT/time", "/PRODUCT/delta_time"]], "VARNAMES": ["/PRODUCT/formaldehyde_tropospheric_vertical_column", "/PRODUCT/qa_value", "/PRODUCT/time_utc", "/PRODUCT/scanline", "/PRODUCT/ground_pixel"], "BOXLONRANGE": [73.0, 136.0], "BOXLATRANGE": [3.0, 54.0], "TIMERANGE": [800414432.0, 800496009.0], "GRIDTYPES": ["SWATH"], "CONVERTFILETYPE": [true]}
    dimensions(sizes): 
    variables(dimensions): 
    groups: PRODUCT, METADATA
In[30]: ns['PRODUCT']
Out[30]: <class 'netCDF4._netCDF4.Group'>
group /PRODUCT:
    dimensions(sizes): time(1), scanline(725), ground_pixel(237)
    variables(dimensions): uint16 time_idx(time), uint16 scanline_idx(scanline), uint16 ground_pixel_idx(ground_pixel), float32 longitude(time,scanline,ground_pixel), float32 latitude(time,scanline,ground_pixel), int32 time(time), int32 delta_time(time,scanline,ground_pixel), float32 formaldehyde_tropospheric_vertical_column(time,scanline,ground_pixel), uint8 qa_value(time,scanline,ground_pixel), <class 'str'> time_utc(time,scanline), int32 scanline(scanline), int32 ground_pixel(ground_pixel)
    groups: SUPPORT_DATA
In[31]: ns['PRODUCT'].variables.keys()
Out[31]: dict_keys(['time_idx', 'scanline_idx', 'ground_pixel_idx', 'longitude', 'latitude', 'time', 'delta_time', 'formaldehyde_tropospheric_vertical_column', 'qa_value', 'time_utc', 'scanline', 'ground_pixel'])

(2) xarray的讀取結果:

xs
Out[34]: 
<xarray.Dataset>
Dimensions:                                    (ground_pixel: 237, scanline: 725, time: 1)
Coordinates:
  * time                                       (time) datetime64[ns] 2018-05-14
  * scanline                                   (scanline) float64 1.507e+03 ....
  * ground_pixel                               (ground_pixel) float64 1.0 ......
Data variables:
    time_idx                                   (time) float32 0.0
    scanline_idx                               (scanline) float32 1.506e+03 ....
    ground_pixel_idx                           (ground_pixel) float32 0.0 ......
    longitude                                  (time, scanline, ground_pixel) float32 ...
    latitude                                   (time, scanline, ground_pixel) float32 ...
    delta_time                                 (time, scanline, ground_pixel) timedelta64[ns] ...
    formaldehyde_tropospheric_vertical_column  (time, scanline, ground_pixel) float32 ...
    qa_value                                   (time, scanline, ground_pixel) float32 ...
    time_utc                                   (time, scanline) object nan .....