目录
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包打开。
总结
以上为个人经验，希望能给大家一个参考，也希望大家多多支持。
© 版权声明
THE END
暂无评论内容