使用xarray更改坐标系以便进行切片操作

问题描述 投票:2回答:1

我是新来的。首先,我非常感谢您的时间和考虑。关于在python中管理2个不同的netcdf文件,我有2个问题。我搜索了很多但不幸的是我找不到解决方案。

1-我有一个netcdf文件,其坐标如下:

time     datetime64[ns] 2016-08-16T22:00:00
* y        (y) int32 220000  ...  620000
* x        (x) int32 20000  ...  720000
 lat      (y, x) float64 dask.array<shape=(401, 701), 
 lon      (y, x) float64 dask.array<shape=(401, 701),

我需要将coords更改为lon / lat,以便我可以根据特定的lon / lat坐标切片区域(通过使用xarray)。但我不知道如何将x和y改为lon lat。我的代码:

import xarray as xr
import matplotlib.pyplot as plt
p = "R_201608.nc"
ds = xr.open_mfdataset(p)
q=ds.RR.sel(time='2016-08-16T21:00:00')

2-与1类似,我有另一个netcdf文件,其坐标如下:

   * X           (X) float32 557600.0 .. 579400.0
   * Y           (Y) float32 5190600 ... 5205400.0
   * time        (time) datetime64[ns] 2007-01I

如何将x和y转换为lon / lat系统,以便我可以在lon / lat系统中绘制它?

编辑与@Ryan相关:1-是。这个文件可以在很大范围内降雨。我想将它切割成更小的区域 - 与q2-相关的文件的相似区域,并使用偏差,RMSE等比较它们。这里是与此文件相关的完整信息:

 <xarray.Dataset>
  Dimensions:                  (time: 2976, x: 701, y: 401)
  Coordinates:
  * time             (time) datetime64[ns] 2016-08-31T23:45:00
  * y          (y) int32 220000 221000  ... 619000 620000
  * x          (x) int32 20000 21000  ... 719000 720000
  lat        (y, x) float64 dask.array<shape=(401, 701),chunksize=(401, 701)>
  lon        (y, x) float64 dask.array<shape=(401, 701), chunksize=(401, 701)

 Data variables:
    RR       (time, y, x) float32 dask.array<shape=(2976, 401, 701),    chunksize=(2976, 401, 701)>
    lambert_conformal_conic  int32 ...

    Conventions:  CF-1.5

与@Ryan相关的编辑:2-这里是关于第二个文件(较小区域)的完整信息:

   <xarray.DataArray 'Precip' (time: 8928, Y: 75, X: 110)>
   dask.array<shape=(8928, 75, 110), dtype=float32, chunksize=(288, 75, 110)>
   Coordinates:

      sensor_height_precip  float32 1.5
      sensor_height_P       float32 1.5
      * X                     (X) float32 557600.0 557800.0 ... 579200.0 579400.0
      * Y                     (Y) float32 5190600.0 5190800.0 ... 5205400.0
      * time                  (time) datetime64[ns]  2007-01-31T23:55:00
   Attributes:
      grid_mapping:         UTM33N
      ancillary_variables:  QFlag_Precip QGrid_Precip
      long_name:            Precipitation Amount
      standard_name:        precipitation_amount
      cell_methods:         time:sum
      units:                mm
python-3.x netcdf coordinate-transformation python-xarray netcdf4
1个回答
1
投票

在问题1)中,不可能将lon和lat转换为维度坐标,因为它们是二维的(都具有维度x,y)。用于切片的尺寸坐标只能是一维的。如果您可以更加具体地了解切片后您想要做什么,我们可以提供更多有关如何进行的建议。是否要选择特定的纬度/经度范围,然后计算一些统计数据(例如均值/方差)?

问题2)看起来你有一个地图投影。如果没有关于投影的更多信息,则无法转换为纬度/经度坐标或在地图上绘图。您的数据集中是否包含有关所使用的地图投影的更多信息?你能发布print(ds)的全部输出吗?


0
投票

我帮你解决了问题。非常感谢。我可以使用PYPROJ将两个数据集的坐标更改为lon / lat,如@Bart所述。从原始坐标和投影坐标创建meshgid是关键点。

from pyproj import Proj
nxv,  nyv = np.meshgrid(nx, ny)       
unausp = Proj('+proj=lcc +lat_1=49 +lat_2=46 +lat_0=47.5   +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +ellps=bessel    +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs ')   
nlons, nlats = unausp(nxv, nyv, inverse=True)                                 
upLon,  upLat = np.meshgrid(nlons,nlats)

由于我想比较具有不同空间分辨率(不同网格大小)的两个降雨数据集,我必须使用xarray插值来升级其中一个:

upnew_lon = np.linspace(w.X[0], w.X[-1], w.dims['X'] // 5) 
upnew_lat = np.linspace(w.Y[0], w.Y[-1], w.dims['Y'] //5) 
uppds = w.interp(Y=upnew_lat, X=upnew_lon)  

据我所知,这种插值是基于线性插值的。我将放大的数据集与原始数据集进行了比较。升级后降雨量平均减少约0.03毫米/天。我只是想知道你认为这种小时降雨的升级方法是否可靠?

© www.soinside.com 2019 - 2024. All rights reserved.