pythonXarray处理设置二维数组作为coordinates方式

Lark ·
更新时间:2024-11-01
· 1366 次阅读

目录

python Xarray处理设置二维数组作为coordinates

Xarray(python)读取​Sentinel-5P(S5P)哨兵数据

使用panoly可视化

使用python里的工具包读取

不足使用xarray读取含Groups的嵌套文件如.nc4时

总结

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 # 创建 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文件的常规包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 ..... 不足使用xarray读取含Groups的嵌套文件如.nc4时

需要先知道其所在的Gropus名称,即需要先用panoly软件或nc4包打开。

总结

以上为个人经验,希望能给大家一个参考,也希望大家多多支持软件开发网。



二维数组 维数 二维 数组

需要 登录 后方可回复, 如果你还没有账号请 注册新账号