Using xarray to change coordinate system in order to Slice operation

Question

I am new here. first on all, I am very thankful for your time and consideration. I have 2 questions regarding to managing 2 different netcdf files in python. I searched a lot but unfortunately I couldn't find a solution.

1- I have a netcdf file which has coordinates like below:

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

I need to change coords to lon/lat in order that I can slice an area based on specific lon/lat coords (by using xarray). But I don't know how to change x and y to lon lat. here my code:

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- Similar to 1, I have another netcdf file which has coordinates like below:

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

How can I convert x and y to lon/lat system in order that I can plot it in lon/lat system?

Edit related to @Ryan : 1- Yes. this file demonestrates rainfall over a large area. I want to cut it into smaller area -similar area of file related to q2- and compare them uusing bias, RMSE, etc. here is full information related to this file:

 <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

edit related to @Ryan :2- And here it is the full information about the second file (Smaller area):

   <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

You can use pyproj for converting the coordinates (your second question). An example for your projection is here. — Bart
– Bart, Commented Mar 17, 2019 at 12:07
Thanks a lot, Bart. I read the example and also searched a lot, but since my file is multidimensional (x,y,t,rain), I couldn't do conversion. — E Ghaemi
– E Ghaemi, Commented Mar 18, 2019 at 16:20
I think I didn't explain the problem well. I have two different netcdf files with different coordination systems and different grid sizes (which I described above). I need to change the coordinate systems so I can regrid one of them and then compare them using some statistical indices (let's say, bias RMSE ...). — E Ghaemi
– E Ghaemi, Commented Mar 19, 2019 at 17:31

Ryan · Accepted Answer · 2019-03-17 01:20:31Z

1

In problem 1), it is not possible to convert lon and lat to dimension coordinates, because they are two-dimensional (both have dimension x, y). Dimension coordinates, used for slicing, can only be one-dimensional. If you can be more specific about what you want to do after slicing, we can provide more suggestions about how to proceed. Do you want to select a particular latitude / longitude range and then calculate some statistics (e.g. mean / variance)?

In problem 2) it looks like you have a map projection. Without more information about the projection, it is impossible to convert to lat / lon coordinates or plot on a map. Is there more information contained in your dataset about the map projection used? Can you post the full output of print(ds)?

answered Mar 17, 2019 at 1:20

Ryan

8067 silver badges13 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

E Ghaemi Over a year ago

Wow, That was quick! Thanks a lot. I added an edition into my question. Thank you very much.

E Ghaemi · Accepted Answer · 2019-06-29 16:39:12Z

I have solved my problem with your help. Thanks a lot. I could change the coords of both data sets to lon/lat using PYPROJ as @Bart mentioned. creating meshgid from original and projected coordinates was the key point.

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)

Since I want to compare two rainfall data sets with different spatial resolution (different grid size), I have to upscale one of them by using xarray interpolation:

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)

AS far as I know, this interpolation is based on linear interpolation. I compared upscaled data set with the original one. The mean of rainfall decreases about 0.03mm/day after upscaling. I just want to know do you think this upscaling method for sub-hourly rainfall is reliable?

Collectives™ on Stack Overflow

Using xarray to change coordinate system in order to Slice operation

2 Answers 2

1 Comment

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Related