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
# 创建 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包打开。
总结
以上为个人经验,希望能给大家一个参考,也希望大家多多支持脚本之家。
相关文章
matplotlib命令与格式之tick坐标轴日期格式(设置日期主副刻度)
这篇文章主要介绍了matplotlib命令与格式之tick坐标轴日期格式(设置日期主副刻度),文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧2019-08-08
jupyter notebook 添加kernel permission denied的操作
这篇文章主要介绍了jupyter notebook 添加kernel permission denied的操作,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧2020-04-04
解决python报错ImportError:urllib3 v2.0 only supports OpenSSL
这篇文章主要介绍了解决python报错ImportError:urllib3 v2.0 only supports OpenSSL 1.1.1+的相关资料,文中通过代码介绍的非常详细,需要的朋友可以参考下2023-12-12
Python+Selenium使用Page Object实现页面自动化测试
这篇文章主要介绍了Python+Selenium使用Page Object实现页面自动化测试,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧2019-07-07


最新评论