ARM-DOE / pyart

The Python-ARM Radar Toolkit. A data model driven interactive toolkit for working with weather radar data.
https://arm-doe.github.io/pyart/
Other
501 stars 264 forks source link

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

Open hanzhe910 opened 3 months ago

hanzhe910 commented 3 months ago

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 commented 3 months ago

@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 commented 3 months ago

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 commented 3 months ago

@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 commented 3 months ago

@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 commented 3 months ago

@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 commented 3 months ago

@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 commented 3 months ago

@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