Использование xarray для изменения системы координат для операции Slice - PullRequest
Новая
       26

Использование xarray для изменения системы координат для операции Slice

3 голосов
/ 17 марта 2019

Я новичок здесь. Прежде всего, я очень благодарен за ваше время и внимание. У меня есть 2 вопроса, касающихся управления 2 различными файлами Netcdf в Python. Я много искал, но, к сожалению, не смог найти решение.

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),

Мне нужно изменить координаты на 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

Ответы [ 2 ]

1 голос
/ 17 марта 2019

В задаче 1) невозможно преобразовать lon и lat в размерные координаты, потому что они двумерные (оба имеют размерность x, y). Координаты размеров, используемые для нарезки, могут быть только одномерными. Если вы можете более конкретно рассказать о том, что вы хотите сделать после нарезки, мы можем предоставить больше советов о том, как действовать дальше. Вы хотите выбрать определенный диапазон широты / долготы, а затем рассчитать статистику (например, среднее значение / дисперсия)?

В задаче 2) похоже, что у вас есть проекция карты. Без дополнительной информации о проекции невозможно преобразовать координаты широты / долготы или построить график на карте. Есть ли в вашем наборе данных дополнительная информация об используемой проекции карты? Можете ли вы опубликовать полный вывод print(ds)?

0 голосов
/ 29 июня 2019

Я решил мою проблему с вашей помощью.Большое спасибо.Я мог бы изменить координаты обоих наборов данных на lon / lat, используя PYPROJ, как упомянуто @Bart.создание сетки из исходных и проектируемых координат было ключевой точкой. 

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 мм / день после повышения.Я просто хочу знать, как вы думаете, надежен ли этот метод масштабирования для ежечасных осадков?

...