Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

An error occurred while reading CAPPI file:ValueError: negative dimensions are not allowed #1585

Open
hanzhe910 opened this issue May 22, 2024 · 7 comments

Comments

@hanzhe910
Copy link

I'm very sorry, I have a very simple question that I hope someone can help me with, I'm trying to open a CAPPI file with pyart, but my code is reporting an error as follows:

ValueError Traceback (most recent call last)
Cell In[11], line 9
3 #pyart.map.grid
5 radar_file_path = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/CAP_V_030_256_20180916110009.CAPSPTC'
----> 9 radar = pyart.io.read(radar_file_path)
12 display = pyart.graph.RadarDisplay(radar)
15 plt.figure(figsize=(10, 6))

File ~/anaconda3/envs/pyart/lib/python3.10/site-packages/pyart/io/auto_read.py:125, in read(filename, use_rsl, **kwargs)
123 return read_rsl(filename, **kwargs)
124 else:
--> 125 return read_sigmet(filename, **kwargs)
126 if filetype == "UF":
127 if use_rsl:

File ~/anaconda3/envs/pyart/lib/python3.10/site-packages/pyart/io/sigmet.py:147, in read_sigmet(filename, field_names, additional_metadata, file_field_names, exclude_fields, include_fields, time_ordered, full_xhdr, noaa_hh_hdr, debug, ignore_xhdr, ignore_sweep_start_ms, **kwargs)
144 full_xhdr = False
146 # read data
--> 147 sigmet_data, sigmet_metadata = sigmetfile.read_data(full_xhdr=full_xhdr)
148 first_data_type = sigmetfile.data_type_names[0]
149 if first_data_type == "XHDR": # don't use XHDR as the first data type

File ~/anaconda3/envs/pyart/lib/python3.10/site-packages/pyart/io/_sigmetfile.pyx:133, in pyart.io._sigmetfile.SigmetFile.read_data()

File ~/anaconda3/envs/pyart/lib/python3.10/site-packages/numpy/ma/core.py:8442, in _convert2ma.call(self, *args, **params)
8440 _extras[p] = params.pop(p)
8441 # Get the result
-> 8442 result = self._func.call(*args, **params).view(MaskedArray)
8443 if "fill_value" in common_params:
8444 result.fill_value = _extras.get("fill_value", None)

ValueError: negative dimensions are not allowed

========================================================================================
Meanwhile, my code is:

import matplotlib.pyplot as plt
import pyart
radar_file_path = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/CAP_V_030_256_20180916110009.CAPSPTC'
radar = pyart.io.read(radar_file_path)
display = pyart.graph.RadarDisplay(radar)
plt.figure(figsize=(10, 6))
sweep = 0
display.plot('reflectivity', sweep=sweep, title='Radar Reflectivity', colorbar_label='dBZ')
plt.show()

I don't know if this is because my file is from an official organization and they identify it as cappi, or maybe it is because I don't have the TRMM Radar Software Library(RSL) installed. But I do not know how to install RSL, this does not seem to be a python package, could someone please help me, thank you very much!!!!

@zssherman
Copy link
Collaborator

@hanzhe910 Hello! Currently there is no CAPPI functionality in Py-ART. However If the file is gridded, you could try pyart.io.read_grid instead for the file and see if that works? If that doesn't we would have to see where the negative dimensions are coming from?

@hanzhe910
Copy link
Author

see

Thank you very much for your reply. I am thrilled to receive a response from the developer of pyart. My CAPPI file is not in netCDF4 format, so it doesn't have a grid, but that's okay. I have communicated with the relevant personnel and am preparing to process the RAW file. Pyart is the best radar data processing package I have ever used. Now, I would like to convert radar data into gridded data. Could you please let me know what formula is used in pyart to convert upper-level winds to surface winds?

@hanzhe910
Copy link
Author

@zssherman
Hello,

I have another question that I would greatly appreciate your help with. Due to the elevation angle coverage of my radar, three files are generated:

PPIVOL_A: 0.1 deg
PPIVOL_B: 0.9, 1.8, 2.7, 3.6, 5.4 deg
PPIVOL_C: 10, 15, 22, 34 deg
I would like to obtain the wind speed at an altitude of 3 km. Therefore, my idea is to combine these three files together. My proposed steps are as follows:

Read the three raw files.
Correct the wind speeds in the three files.
Grid the three files to the 3 km altitude.
Merge the three gridded files together.
Could you please let me know if this approach is feasible? I have noticed that all these steps can be accomplished using pyart. However, I have another question: when I use pyart.map.grid_from_radars() to grid the three files, does pyart take into account the tilt angles of the wind speeds?

Thank you very much for your assistance!!!!!!

@zssherman
Copy link
Collaborator

@hanzhe910 In terms of converting upper level winds to surface winds, I'm not entirely sure. Someone else might have an idea on that. For joining radars, we do have a pyart.utils.join_radar that will allow you to join radars together. You could join them by the order of the sweeps so ppi_volA with B then those two with C i think would work. You could also use the pyart.correct.region_delias to correct velocities. The pyart horizontal wind profile class could help possibly as well if i understand the problem correctly. Or could use the vad_browning pyart.retrieve function and specify the height you want to obtain a VAD ?

@hanzhe910
Copy link
Author

@zssherman
Thank you very, very much for your detailed and helpful response. I am extremely grateful for your assistance and the valuable insights you provided. Your suggestions have been incredibly beneficial, and I believe I can use the pyart.retrieve.vad_michelson function to perform the VAD analysis and address my issue.

However, I have one more question. When reading the radar wind data, I understand that it refers to the radial wind, which forms an angle with the ground and is not the horizontal wind. How can I obtain the horizontal wind from this data? Does Py-ART provide functionality to achieve this, or can the projection feature in pyart.map.map_to_grid be used for this purpose?

Thank you so much once again for your invaluable help and prompt response.

@hanzhe910
Copy link
Author

hanzhe910 commented May 23, 2024

@zssherman

I have successfully merged radar data and applied the dealias_region_based function, followed by gridding using pyart.map.grid_from_radars. Here is my code:

==============================================================================
import os
import matplotlib.pyplot as plt
import pyart
import numpy as np
import xarray as xr

radar_file_path1 = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar1'
radar_file_path2 = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar2'
radar_file_path3 = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar3'

radar1 = pyart.io.read(radar_file_path1)
radar2 = pyart.io.read(radar_file_path2)
radar3 = pyart.io.read(radar_file_path3)

radarcomb = pyart.util.join_radar(radar1, radar2)
radarcombfinal = pyart.util.join_radar(radarcomb, radar3)

output_dir = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar_plots'
os.makedirs(output_dir, exist_ok=True)

def plot_radar_sweeps(radar, radar_name, output_dir):
display = pyart.graph.RadarDisplay(radar)
num_sweeps = radar.nsweeps
for sweep in range(num_sweeps):
plt.figure(figsize=(10, 6))
display.plot('velocity',vmin=-20,vmax=60, sweep=sweep, title=f'{radar_name} Sweep {sweep + 1}', colorbar_label='m/s')
plt.savefig(os.path.join(output_dir, f'{radar_name}sweep{sweep + 1}.png'))
plt.close()

plot_radar_sweeps(radarcombfinal, 'Radar', output_dir)
print(f'All plots are saved in {output_dir}')

dealias_output_dir = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar_plots_dealias'
os.makedirs(dealias_output_dir, exist_ok=True)

dealiased_velocity = pyart.correct.dealias_region_based(radarcombfinal)

radarcombfinal.add_field('velocity', dealiased_velocity, replace_existing=True)

plot_radar_sweeps(radarcombfinal, 'Radar_Dealiased', dealias_output_dir)
print(f'All dealiased plots are saved in {dealias_output_dir}')

z_grid_limits = (3000,3150)
y_grid_limits = (-150000.,150000.)
x_grid_limits = (-150000.,150000.)

grid_resolution = 150

def compute_number_of_points(extent, resolution):
return int((extent[1] - extent[0])/resolution)

z_grid_points = compute_number_of_points(z_grid_limits, grid_resolution)
x_grid_points = compute_number_of_points(x_grid_limits, grid_resolution)
y_grid_points = compute_number_of_points(y_grid_limits, grid_resolution)

print(z_grid_points,
y_grid_points,
x_grid_points)

grid = pyart.map.grid_from_radars(radarcombfinal,
grid_shape=(z_grid_points,
y_grid_points,
x_grid_points),
grid_limits=(z_grid_limits,
y_grid_limits,
x_grid_limits))
#weighting_function='Barnes2',
#gridding_algo='map_to_grid')

display = pyart.graph.GridMapDisplay(grid)
display.plot_grid('velocity',
level=0,
vmin=-20,
vmax=60,
cmap='pyart_HomeyerRainbow')

gridxarray = grid.to_xarray()

xarray_output_dir = '/Users/cuihanzhe/Desktop/ai_downscaling/radardata/radar_plots_xarray'
os.makedirs(xarray_output_dir, exist_ok=True)

plt.figure(figsize=(10, 6))

gridxarray.isel(z=0).velocity.plot(cmap='pyart_HomeyerRainbow',
vmin=-20,
vmax=60);

plt.title(f'Wind Speed at 3 km')
plt.savefig(os.path.join(xarray_output_dir, f'wind_speed_3km.png'))
plt.close()

======================================================================================
However, I encountered a significant issue with my gridding results. In the resulting image after gridding, there are concentric circles at the center where the data is missing, showing as white concentric circles.
wind_speed_3km
I have verified that my radar data does not have this issue because the radar images processed with dealias_region_based are perfectly normal without any concentric circles.
Radar_Dealiased_sweep_1
Radar_Dealiased_sweep_5
Radar_Dealiased_sweep_10

Could you please help me understand why this might be happening? I would greatly appreciate any assistance you can provide.

Thank you very much for your help!

@zssherman
Copy link
Collaborator

@hanzhe910 Hmm I would maybe try tweaking the ROI (but might not be ideal) of the gridding etc. But if I had to take a guess it might be something along the lines of this or the link could be helpful in pinpointing why that's the case:
https://groups.google.com/g/pyart-users/c/9G5T9Q9wrPA/m/VIVpIqL1AgAJ

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants