openradar / xradar

A tool to work in weather radar data in xarray
https://docs.openradarscience.org/projects/xradar
MIT License
85 stars 17 forks source link

`pyproj.exceptions.CRSError` when georeferencing `cfradial1` file #146

Closed aladinor closed 8 months ago

aladinor commented 8 months ago

This question is related to #145 . After open a Cfradial file I tried to apply georeference as follows

import xradar as xd
import numpy as np
from open_radar_data import DATASETS

# filepath = DATASETS.fetch('cfrad.20211011_201733.023_to_20211011_201745.299_DOW8_RHI.nc')
filepath = DATASETS.fetch('cfrad.20211011_201711.345_to_20211011_201732.860_DOW8_PPI.nc')
dt = xd.io.open_cfradial1_datatree(filepath)
ds = dt['sweep_0'].ds.copy()
ds = ds.xradar.georeference()

Then, I got a pyproj.exceptions.CRSError

  File "pyproj/_crs.pyx", line 2378, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection: +proj=aeqd +datum=WGS84 +lon_0=[-88.33178711 -88.33178711 -88.33178711 -88.33178711 -88.33178711
 -88.33178711 -88.33178711 -88.33178711 -88.33178711 -88.33178711
 -88.33178711 -88.33178711 -88.33178711 -88.33178711 -88.33178711

Digging a little bit, this happens because this cfradial file reports more than 1 value in the longitude, latitude and elevation dimensions/variables. Therefore, when creating the projparams, the pyproj.exceptions.CRSError raise up. To solve this issue, I found out that I can take the mean/median value for those variables, and everything works perfectly.

ds = dt['sweep_0'].ds.copy()
ds['longitude'] = ds['longitude'].mean('time').values
ds['latitude'] = ds['latitude'].mean('time').values
ds['azimuth'] = np.median(ds['azimuth'].values)
ds['altitude'] = np.median(ds['altitude'].values)
ds = ds.xradar.georeference()

If we look at the ds

<xarray.Dataset>
Dimensions:                    (time: 820, range: 1984, azimuth: 820)
Coordinates:
  * time                       (time) datetime64[ns] 2021-10-11T20:17:11.3450...
  * range                      (range) float32 37.47 112.4 ... 1.487e+05
    elevation                  (azimuth) float32 0.4779 0.4779 ... 0.6042 0.6042
    latitude                   float64 40.01
    longitude                  float64 -88.33
    altitude                   float64 212.0
  * azimuth                    (azimuth) float32 33.84 33.84 ... 326.1 326.2
    crs_wkt                    int64 0
    x                          (azimuth, range) float32 20.87 ... -8.276e+04
    y                          (azimuth, range) float32 31.12 ... 1.234e+05
    z                          (azimuth, range) float32 212.0 213.0 ... 3.08e+03
Data variables: (12/31)
    sweep_number               float64 ...
    sweep_mode                 <U6 'sector'
    prt_mode                   |S32 ...
    follow_mode                |S32 ...
    sweep_fixed_angle          float32 ...
    ray_start_range            (azimuth) float32 ...
    ...                         ...
    DBMHC                      (azimuth, range) float32 ...
    DBZHC                      (azimuth, range) float32 ...
    VEL                        (azimuth, range) float32 ...
    VS1                        (azimuth, range) float32 ...
    VL1                        (azimuth, range) float32 ...
    WIDTH                      (azimuth, range) float32 ...

Should we care about this in the Cfradial accessor stage? or, where do you think we could take care of this?

Cheers,

Alfonso.

kmuehlbauer commented 8 months ago

@aladinor Can you show the output of the original ds, before adapting the coords?

The question is why would there be a change in coordinates? Does the DOW move while scanning?

syedhamidali commented 8 months ago

@kmuehlbauer I think it is moving a little from the below-given info, and so is its name suggesting, "Doppler on Wheels (DOW)"

np.unique(ds.longitude.values), np.unique(ds.latitude.values), np.unique(ds.altitude.values)
(array([-88.33179474, -88.33178711,          nan]),
 array([40.0148201 , 40.01482391,         nan]),
 array([210.99999547, 211.9999975 ,          nan]))
syedhamidali commented 8 months ago

@aladinor i overlooked your issue, and didn't notice you're also providing solution to it.

aladinor commented 8 months ago

Hi @kmuehlbauer and @syedhamidali. I didn't know that DOW8 data records lat, lon, and elevation as a function of time. This may be produced as it is a mobile platform. I thought it was an I/O error, but it isn't. I will close this issue since it is not properly a bug. Thanks once more for your time and help.

Cheers,

Alfonso