radar_data.py:
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.transform import from_origin
from rasterio.warp import calculate_default_transform, reproject, Resampling
import geopandas as gpd
from shapely.geometry import Point, Polygon
from pyproj import Proj, Transformer
import struct
from tqdm import tqdm
import os
from lxml import etree
import io
import sys
import logging

# Configure logging
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s',
    handlers=[
        logging.StreamHandler(sys.stdout),
        logging.FileHandler('radar_data_processing.log', mode='w')
    ]
)

# Define constants and paths
RADAR_DATA_PATH = 'input/radar/dBZ.vol'  # Path to the radar data file
MASOVIAN_BOUNDARY_PATH = 'input/radar/masovian_boundary.geojson'  # Path to Masovian boundary GeoJSON or Shapefile
OUTPUT_PLOT_PATH = 'output/processed_radar_data.png'  # Path to save the output plot

# Radar specifications based on provided XML
AZIMUTH_RESOLUTION = 1.0  # Degrees between each azimuth
RANGE_RESOLUTION = 0.5  # As per <rangestep>0.5</rangestep> km

# Radar location from <sensorinfo>
RADAR_LATITUDE = 52.405250  # From <lat>52.405250</lat>
RADAR_LONGITUDE = 20.961110  # From <lon>20.961110</lon>
RADAR_ALTITUDE = 122.3  # From <alt>122.300000</alt> meters

# Coordinate Systems
SOURCE_CRS = 'EPSG:4326'  # Geographic coordinates
TARGET_CRS = 'EPSG:2180'  # Polish coordinate system

def parse_vol_file(file_path):
    """
    Parses the .vol file to extract metadata about each radar slice.
    
    Parameters:
    - file_path: Path to the .vol file.
    
    Returns:
    - slices_info: List of dictionaries containing slice metadata.
    - sensor_info: Dictionary containing radar sensor information.
    - binary_data_start: Byte index where binary data starts.
    """
    logging.info("Parsing the .vol file...")
    slices_info = []
    sensor_info = {}
    
    with open(file_path, 'rb') as f:
        # Read the entire file
        content = f.read()
    
    # Locate the end of the XML content
    end_volume_tag = b'</volume>'
    end_volume_index = content.find(end_volume_tag)
    if end_volume_index == -1:
        logging.error("</volume> tag not found in the file.")
        sys.exit(1)
    
    xml_content = content[:end_volume_index + len(end_volume_tag)]
    binary_content = content[end_volume_index + len(end_volume_tag):]
    
    # Parse XML
    try:
        root = etree.fromstring(xml_content)
    except etree.XMLSyntaxError as e:
        logging.error(f"Failed to parse XML content: {e}")
        sys.exit(1)
    
    # Extract sensor info
    sensor = root.find('.//sensorinfo')
    if sensor is not None:
        sensor_info['id'] = sensor.findtext('id')
        sensor_info['type'] = sensor.findtext('type')
        sensor_info['name'] = sensor.findtext('name')
        try:
            sensor_info['lon'] = float(sensor.findtext('lon'))
            sensor_info['lat'] = float(sensor.findtext('lat'))
            sensor_info['alt'] = float(sensor.findtext('alt'))
            sensor_info['wavelen'] = float(sensor.findtext('wavelen'))
            sensor_info['beamwidth'] = float(sensor.findtext('beamwidth'))
        except (TypeError, ValueError) as e:
            logging.error(f"Invalid sensor information: {e}")
            sys.exit(1)
    else:
        logging.error("<sensorinfo> not found in the XML.")
        sys.exit(1)
    
    # Extract slice information
    slices = root.findall('.//slice')
    for slice_elem in slices:
        slice_info = {}
        slice_info['refid'] = slice_elem.get('refid')
        # Find rawdata element for dBZ
        rawdata = slice_elem.find('.//rawdata[@type="dBZ"]')
        if rawdata is not None:
            try:
                slice_info['rays'] = int(rawdata.get('rays'))
                slice_info['bins'] = int(rawdata.get('bins'))
                slice_info['type'] = rawdata.get('type')
                slice_info['max'] = float(rawdata.get('max'))
                slice_info['min'] = float(rawdata.get('min'))
                slices_info.append(slice_info)
                logging.info(f"Slice refid {slice_info['refid']}: rays {slice_info['rays']}, bins {slice_info['bins']}.")
            except (TypeError, ValueError) as e:
                logging.warning(f"Invalid slice information in refid {slice_info.get('refid')}: {e}. Skipping this slice.")
        else:
            logging.warning(f"No dBZ rawdata found for slice refid {slice_elem.get('refid')}. Skipping this slice.")
    
    # Determine where binary data starts
    binary_data_start = end_volume_index + len(end_volume_tag)
    
    logging.info(f"Total slices with dBZ data: {len(slices_info)}.\n")
    
    return slices_info, sensor_info, binary_data_start, binary_content

def extract_dBZ_data(slices_info, binary_content, binary_data_start):
    """
    Extracts dBZ data from the binary content based on slice metadata.
    
    Parameters:
    - slices_info: List of dictionaries containing slice metadata.
    - binary_content: Bytes object containing all binary data after </volume>.
    - binary_data_start: Byte index where binary data starts.
    
    Returns:
    - dBZ_combined: 2D numpy array combining all slices (azimuths x range_bins).
    """
    logging.info("Extracting dBZ data from binary content...")
    dBZ_list = []
    offset = 0  # Current byte offset within binary_content
    
    for idx, slice_info in enumerate(slices_info):
        rays = slice_info['rays']
        bins = slice_info['bins']
        dtype = np.float32  # float32 per dBZ value
        byte_size = 4  # bytes per float32
        expected_bytes = rays * bins * byte_size
        
        # Check if enough data remains
        if offset + expected_bytes > len(binary_content):
            logging.warning(f"Slice {idx + 1} refid {slice_info['refid']}: Insufficient binary data. Expected {expected_bytes} bytes, got {len(binary_content) - offset} bytes. Adjusting number of bins.")
            # Adjust the number of bins based on remaining data
            available_bytes = len(binary_content) - offset
            valid_bins = available_bytes // (rays * byte_size)
            if valid_bins <= 0:
                logging.error(f"Slice {idx + 1} refid {slice_info['refid']}: No data available. Skipping this slice.")
                continue
            slice_info['bins'] = valid_bins
            bins = valid_bins
            expected_bytes = rays * bins * byte_size
        
        # Extract binary data for this slice
        slice_binary_data = binary_content[offset:offset + expected_bytes]
        offset += expected_bytes  # Move to the next slice
        
        try:
            # Convert binary data to numpy array
            dBZ = np.frombuffer(slice_binary_data, dtype=dtype).reshape((rays, bins))
            dBZ_list.append(dBZ)
            logging.info(f"Slice {idx + 1} refid {slice_info['refid']}: Extracted dBZ data with shape {dBZ.shape}.")
        except Exception as e:
            logging.error(f"Slice {idx + 1} refid {slice_info['refid']}: Failed to parse dBZ data. Error: {e}. Skipping this slice.")
            continue
    
    if not dBZ_list:
        logging.error("No dBZ data extracted. Exiting.")
        sys.exit(1)
    
    # Ensure all slices have the same number of bins by padding with np.nan
    max_bins = max(dBZ.shape[1] for dBZ in dBZ_list)
    logging.info(f"Maximum number of bins across slices: {max_bins}")
    
    padded_dBZ_list = []
    for idx, dBZ in enumerate(dBZ_list):
        current_bins = dBZ.shape[1]
        if current_bins < max_bins:
            padding = ((0, 0), (0, max_bins - current_bins))
            dBZ_padded = np.pad(dBZ, padding, mode='constant', constant_values=np.nan)
            logging.info(f"Padded slice {idx + 1} from {current_bins} to {max_bins} bins.")
        else:
            dBZ_padded = dBZ
        padded_dBZ_list.append(dBZ_padded)
    
    # Stack all slices vertically
    try:
        dBZ_combined = np.vstack(padded_dBZ_list)
        logging.info(f"Combined dBZ data shape after padding and stacking: {dBZ_combined.shape}\n")
    except Exception as e:
        logging.error(f"Failed to stack dBZ data: {e}")
        sys.exit(1)
    
    return dBZ_combined

def polar_to_cartesian(dBZ, radar_lat, radar_lon, az_res, rng_res):
    """
    Converts polar radar data to Cartesian coordinates.
    
    Parameters:
    - dBZ: 2D numpy array of dBZ values (azimuths x ranges).
    - radar_lat: Latitude of the radar.
    - radar_lon: Longitude of the radar.
    - az_res: Azimuth resolution in degrees.
    - rng_res: Range resolution in km.
    
    Returns:
    - xs: 1D numpy array of x coordinates in meters.
    - ys: 1D numpy array of y coordinates in meters.
    - dBZ_1D: 1D numpy array of dBZ values.
    """
    logging.info("Converting polar coordinates to Cartesian coordinates...")
    num_az, num_rng = dBZ.shape
    azimuths = np.deg2rad(np.arange(0, 360, az_res))
    ranges = np.arange(rng_res, rng_res * (num_rng + 1), rng_res) * 1000  # Convert km to meters
    
    # Ensure azimuths match the number of rays
    if len(azimuths) != num_az:
        logging.warning(f"Number of azimuths ({len(azimuths)}) does not match expected ({num_az}). Adjusting.")
        azimuths = np.linspace(0, 2 * np.pi, num=num_az, endpoint=False)
    
    # Adjust ranges if num_rng does not match expected
    if len(ranges) != num_rng:
        logging.warning(f"Number of range bins ({len(ranges)}) does not match expected ({num_rng}). Adjusting.")
        ranges = np.linspace(rng_res, rng_res * num_rng, num=num_rng) * 1000  # Ensure in meters
    
    # Create meshgrid
    az, r = np.meshgrid(azimuths, ranges, indexing='ij')
    
    # Convert to Cartesian coordinates
    x = r * np.sin(az)
    y = r * np.cos(az)
    
    # Flatten the arrays
    xs = x.flatten()
    ys = y.flatten()
    dBZ_1D = dBZ.flatten()
    
    logging.info("Coordinate conversion completed.\n")
    return xs, ys, dBZ_1D

def create_geodataframe(xs, ys, dBZ):
    """
    Creates a GeoDataFrame from Cartesian coordinates and dBZ values.
    
    Parameters:
    - xs: 1D numpy array of x coordinates in meters.
    - ys: 1D numpy array of y coordinates in meters.
    - dBZ: 1D numpy array of dBZ values.
    
    Returns:
    - gdf: GeoDataFrame with geometry in geographic coordinates.
    """
    logging.info("Creating GeoDataFrame...")
    # Define local projection (using Transverse Mercator centered at radar location)
    proj_local = Proj(proj='tmerc', lat_0=RADAR_LATITUDE, lon_0=RADAR_LONGITUDE, ellps='WGS84')
    transformer = Transformer.from_proj(proj_local, 'epsg:4326', always_xy=True)
    
    # Transform to geographic coordinates
    lon, lat = transformer.transform(xs, ys)
    
    # Create GeoDataFrame
    geometry = [Point(xy) for xy in zip(lon, lat)]
    gdf = gpd.GeoDataFrame({'dBZ': dBZ}, geometry=geometry, crs='EPSG:4326')
    
    logging.info("GeoDataFrame creation completed.\n")
    return gdf

def reproject_gdf(gdf, target_crs='EPSG:2180'):
    """
    Reprojects the GeoDataFrame to the target CRS.
    
    Parameters:
    - gdf: Input GeoDataFrame.
    - target_crs: Target Coordinate Reference System.
    
    Returns:
    - gdf_reprojected: Reprojected GeoDataFrame.
    """
    logging.info(f"Reprojecting GeoDataFrame to {target_crs}...")
    try:
        gdf_reprojected = gdf.to_crs(target_crs)
        logging.info("Reprojection completed.\n")
    except Exception as e:
        logging.error(f"Failed to reproject GeoDataFrame: {e}")
        sys.exit(1)
    return gdf_reprojected

def clip_to_masovian(gdf, masovian_boundary):
    """
    Clips the GeoDataFrame to the Masovian boundary.
    
    Parameters:
    - gdf: Input GeoDataFrame.
    - masovian_boundary: GeoDataFrame or Polygon of Masovian boundary.
    
    Returns:
    - gdf_clipped: Clipped GeoDataFrame.
    """
    logging.info("Clipping radar data to Masovian boundary...")
    try:
        # Spatial join or intersection
        gdf_clipped = gdf[gdf.within(masovian_boundary.unary_union)]
        logging.info(f"Clipped data contains {len(gdf_clipped)} points.\n")
    except Exception as e:
        logging.error(f"Failed to clip data to Masovian boundary: {e}")
        sys.exit(1)
    return gdf_clipped

def visualize_radar_data(gdf, masovian_boundary, output_path):
    """
    Visualizes the radar data.
    
    Parameters:
    - gdf: GeoDataFrame containing clipped radar data.
    - masovian_boundary: GeoDataFrame or Polygon of Masovian boundary.
    - output_path: Path to save the plot.
    
    Returns:
    - None
    """
    logging.info("Visualizing radar data...")
    try:
        fig, ax = plt.subplots(figsize=(12, 12))
        
        # Plot Masovian boundary
        masovian_boundary.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=1)
        
        # Plot radar dBZ values using scatter plot
        # Filter out invalid dBZ values (e.g., NaN or extreme values)
        valid_gdf = gdf[(gdf['dBZ'] >= 0) & (gdf['dBZ'] <= 80)]
        
        scatter = valid_gdf.plot(
            ax=ax,
            column='dBZ',
            cmap='jet',
            markersize=1,
            alpha=0.6,
            legend=True,
            legend_kwds={'label': "dBZ Value", 'shrink': 0.5}
        )
        
        ax.set_title("Processed Radar dBZ Data within Masovian Region", fontsize=16)
        ax.set_xlabel("Easting (meters)", fontsize=12)
        ax.set_ylabel("Northing (meters)", fontsize=12)
        ax.set_aspect('equal')
    
        # Remove axis ticks for cleaner look
        ax.set_xticks([])
        ax.set_yticks([])
        
        # Add colorbar
        sm = plt.cm.ScalarMappable(cmap='jet', norm=plt.Normalize(vmin=valid_gdf['dBZ'].min(), vmax=valid_gdf['dBZ'].max()))
        sm._A = []  # dummy array for the scalar mappable
        cbar = fig.colorbar(sm, ax=ax, fraction=0.036, pad=0.04)
        cbar.set_label('dBZ', fontsize=12)
    
        # Save the plot
        plt.savefig(output_path, dpi=300, bbox_inches='tight')
        plt.close()
        logging.info(f"Radar data visualization saved to {output_path}.\n")
    except Exception as e:
        logging.error(f"Failed to visualize radar data: {e}")
        sys.exit(1)

def main():
    # Step 1: Parse the .vol file
    if not os.path.exists(RADAR_DATA_PATH):
        logging.error(f"Radar data file not found at {RADAR_DATA_PATH}. Please check the path.")
        sys.exit(1)
    
    slices_info, sensor_info, binary_data_start, binary_content = parse_vol_file(RADAR_DATA_PATH)
    
    # Step 2: Extract dBZ data from binary content
    dBZ_combined = extract_dBZ_data(slices_info, binary_content, binary_data_start)
    
    # Step 3: Convert polar to Cartesian coordinates
    xs, ys, dBZ_1D = polar_to_cartesian(
        dBZ_combined, 
        RADAR_LATITUDE, 
        RADAR_LONGITUDE, 
        AZIMUTH_RESOLUTION, 
        RANGE_RESOLUTION
    )
    
    # Step 4: Create GeoDataFrame
    gdf_geo = create_geodataframe(xs, ys, dBZ_1D)
    
    # Step 5: Reproject to EPSG:2180
    gdf_proj = reproject_gdf(gdf_geo, TARGET_CRS)
    
    # Step 6: Load Masovian Boundary
    if not os.path.exists(MASOVIAN_BOUNDARY_PATH):
        logging.error(f"Masovian boundary file not found at {MASOVIAN_BOUNDARY_PATH}. Please check the path.")
        sys.exit(1)
    try:
        masovian_boundary = gpd.read_file(MASOVIAN_BOUNDARY_PATH)
        logging.info("Loaded Masovian boundary successfully.\n")
    except Exception as e:
        logging.error(f"Failed to load Masovian boundary: {e}")
        sys.exit(1)
    
    # Ensure Masovian boundary is in the same CRS as target CRS
    if masovian_boundary.crs != TARGET_CRS:
        try:
            masovian_boundary = masovian_boundary.to_crs(TARGET_CRS)
            logging.info(f"Reprojected Masovian boundary to {TARGET_CRS}.\n")
        except Exception as e:
            logging.error(f"Failed to reproject Masovian boundary: {e}")
            sys.exit(1)
    
    # Step 7: Clip Radar Data to Masovian Boundary
    gdf_clipped = clip_to_masovian(gdf_proj, masovian_boundary)
    
    # Step 8: Visualize and Plot
    visualize_radar_data(gdf_clipped, masovian_boundary, OUTPUT_PLOT_PATH)
    
    logging.info("Radar data processing and visualization completed successfully.")

if __name__ == "__main__":
    main()
dBZ.vol:
<volume type="vol" version="5.63.6" datetime="2024-10-05T18:49:14Z" owner="">
   <scan date="2024-10-05" datetimehighaccuracy="2024-10-05T18:48:54.193" time="18:48:54" datetimeendofslice="false" name="LEG_125_ZVW.vol">
      <unitid>SI</unitid>
      <advancedchanged>0</advancedchanged>
      <detailedchanged>1</detailedchanged>
      <scantime>249</scantime>
      <starttime>2024-10-05T18:46:32</starttime>
      <endtime>2024-10-05T18:49:14</endtime>
      <nominalstarttime>2024-10-05T18:46:30</nominalstarttime>
      <nominalendtime>2024-10-05T18:51:30</nominalendtime>
      <pargroup refid="sdfbase">
         <multitrip>Off</multitrip>
         <numtrip>2</numtrip>
         <maxsuppressedtrip>3</maxsuppressedtrip>
         <triprecovery>Off</triprecovery>
         <multitripprfmode>Off</multitripprfmode>
         <phase_coding>2</phase_coding>
         <phase_pattern>0</phase_pattern>
         <stoprange>250</stoprange>
         <numele>8</numele>
         <scanstrategy>Optimised</scanstrategy>
         <firstele>0.5</firstele>
         <lastele>40</lastele>
         <antdirection>0</antdirection>
         <sectorscan>Off</sectorscan>
         <gdrx_dp_proc_mode>GdrxDpModeHV_HV</gdrx_dp_proc_mode>
         <highprf>550.00</highprf>
         <stagger>None</stagger>
         <lowprf>0.00</lowprf>
         <midprf>0.00</midprf>
         <nyquistextensionfactor>2.05</nyquistextensionfactor>
         <nyquistextensionmode>FreePRF</nyquistextensionmode>
         <start_range>0</start_range>
         <rangestep>0.5</rangestep>
         <rangesamp>1</rangesamp>
         <cmd>Off</cmd>
         <antspeed>24</antspeed>
         <anglestep>1</anglestep>
         <timesamp>23</timesamp>
         <dynpw_matchedrangeresolution>None</dynpw_matchedrangeresolution>
         <pw_index>0</pw_index>
         <cf_doppler>Off</cf_doppler>
         <fftfilter>On</fftfilter>
         <filtermode>FFT</filtermode>
         <cf_spatial>On</cf_spatial>
         <speckle>Off</speckle>
         <gdrx_clutter_flag_filter>Off</gdrx_clutter_flag_filter>
         <cf_dectree>Off</cf_dectree>
         <cf_gip>Off</cf_gip>
         <matchedfiltertype>NoAlgType</matchedfiltertype>
         <filterwidth>6</filterwidth>
         <cf_gip_width>0.5000</cf_gip_width>
         <cf_gip_mode_const>On</cf_gip_mode_const>
         <gdrx_cluttermap>GdrxCmOff</gdrx_cluttermap>
         <int_filter_mode>1</int_filter_mode>
         <int_filter_corr_mode>Interpolate</int_filter_corr_mode>
         <int_filter_hconf>18</int_filter_hconf>
         <int_filter_vconf>18</int_filter_vconf>
         <int_filter_phaseint>On</int_filter_phaseint>
         <fixselect>AngleStep</fixselect>
         <superresolution>false</superresolution>
         <csr>30</csr>
         <sqi>0.3</sqi>
         <zsqi>0.5</zsqi>
         <log>2</log>
         <extrarclfile></extrarclfile>
         <comment></comment>
         <dualprfmode>SdfDPrfModeAda</dualprfmode>
         <prfcorrthr>2</prfcorrthr>
         <prfcorron>On</prfcorron>
         <prfcorrmode>SdfDPrfCorrMedian</prfcorrmode>
         <filterdepth>SdfCfDepthDefault</filterdepth>
         <fft_weight>GdrxWfHamming</fft_weight>
         <spatial_modus>SdfSpatialFilterAvg</spatial_modus>
         <spatial_min>3</spatial_min>
         <spatial_max>5</spatial_max>
         <gdrx_tripnotchwidth>0.5</gdrx_tripnotchwidth>
         <gdrx_tripccorthres>20</gdrx_tripccorthres>
         <gdrx_tripsnrthres>3</gdrx_tripsnrthres>
         <gdrx_tripsqithres>0.4</gdrx_tripsqithres>
         <gdrx_tripdftwin>GdrxWfRectangular</gdrx_tripdftwin>
         <datatypes>CCOR,CCORv,dBZ,dBuZ,dBZv,dBuZv,FLAGS,uKDP,KDP,uPhiDP,PhiDP,RhoHV,SNR,SNRu,SNRv,SQI,SQIv,V,W,ZDR</datatypes>
         <masterdatatypes>CCOR,dBZ,FLAGS,KDP,PhiDP,RhoHV,SNR,SQI,V,W,ZDR</masterdatatypes>
         <datatypes_format>dBZ=0,dBZv=0,dBuZ=0,dBuZv=0,V=0,Vv=0,Vu=0,Vvu=0,W=0,Wv=0,Wu=0,Wvu=0,ZDR=0,ZDRu=0,uPhiDP=1,uPhiDPu=1,PhiDP=1,KDP=0,uKDP=0,uKDPu=0,RhoHV=0,RhoHVu=0,CCOR=0,CCORv=0,SNR=1,SNRv=0,SNRu=0,SNRvu=0,SQI=0,SQIu=0,SQIv=0,SQIvu=0,FLAGS=-1</datatypes_format>
         <posanglelist>23.800000,15.000000,9.100000,5.400000,3.500000,2.500000,1.500000,0.500000,3.500000,2.500000,1.500000,0.500000</posanglelist>
         <dynpw>0.50</dynpw>
         <dynpw_constantdc_enabled>Off</dynpw_constantdc_enabled>
         <dynpw_constantdc>0.00</dynpw_constantdc>
         <useonlinenoisepower>Off</useonlinenoisepower>
      </pargroup>
      <slice refid="0">
         <dynz max="95.5" min="-31.5"/>
         <dynv max="47.79" min="-47.79"/>
         <dynw max="15.8678" min="0.0622266"/>
         <dynzdr max="12" min="-8"/>
         <dynldr max="10" min="-50"/>
         <dynkdp max="10" min="-10"/>
         <dynsnr max="127" min="0"/>
         <posangle>0.5</posangle>
         <antspeed>18</antspeed>
         <timesamp>29</timesamp>
         <calibslottypes>DATA</calibslottypes>
         <cf_gip_width>0.4</cf_gip_width>
         <int_filter_mode>3</int_filter_mode>
         <int_filter_hconf>10</int_filter_hconf>
         <int_filter_vconf>10</int_filter_vconf>
         <prtstaggering>0</prtstaggering>
         <prtprf>None</prtprf>
         <threshold>dBZ:622602,dBZv:622602,dBuZ:1,dBuZv:1,V:524294,Vv:524294,Vu:1,Vvu:1,W:524294,Wv:524294,Wu:1,Wvu:1,ZDR:622602,ZDRu:1,uPhiDP:2,uPhiDPu:1,PhiDP:524290,KDP:524290,uKDP:2,uKDPu:1,RhoHV:622594,RhoHVu:1</threshold>
         <numtrip>2</numtrip>
         <maxsuppressedtrip>3</maxsuppressedtrip>
         <stoprange>125</stoprange>
         <start_range>0</start_range>
         <rangestep>0.5</rangestep>
         <rangesamp>4</rangesamp>
         <sectorscan>Off</sectorscan>
         <cf_doppler>Off</cf_doppler>
         <fft_weight>GdrxWfHamming</fft_weight>
         <spatial_modus>SdfSpatialFilterAvg</spatial_modus>
         <spatial_min>3</spatial_min>
         <spatial_max>5</spatial_max>
         <startangle>0.0</startangle>
         <stopangle>0.0</stopangle>
         <highprf>1200</highprf>
         <stagger>4/3</stagger>
         <prfcorrmode>SdfDPrfCorrMedian</prfcorrmode>
         <prfcorron>On</prfcorron>
         <prfcorrthr>2</prfcorrthr>
         <dualprfmode>SdfDPrfModeAda</dualprfmode>
         <lowprf>900.00</lowprf>
         <fftfilter>Off</fftfilter>
         <filterwidth>4</filterwidth>
         <filterdepth>SdfCfDepthDefault</filterdepth>
         <statisticalfilter>Off</statisticalfilter>
         <cf_spatial>Off</cf_spatial>
         <superresolution>Off</superresolution>
         <antdirection>0</antdirection>
         <fixselect>AngleStep</fixselect>
         <anglestep>1</anglestep>
         <asyncoffset>0</asyncoffset>
         <anglesync>On</anglesync>
         <csr>30</csr>
         <sqi>0.1</sqi>
         <zsqi>0.1</zsqi>
         <log>2</log>
         <speckle>On</speckle>
         <extrarclfile></extrarclfile>
         <phase_coding>2</phase_coding>
         <gdrx_cluttermap>GdrxCmOff</gdrx_cluttermap>
         <gdrx_tripnotchwidth>0.5</gdrx_tripnotchwidth>
         <gdrx_tripccorthres>20</gdrx_tripccorthres>
         <gdrx_tripsnrthres>3</gdrx_tripsnrthres>
         <gdrx_tripsqithres>0.4</gdrx_tripsqithres>
         <gdrx_tripdftwin>GdrxWfRectangular</gdrx_tripdftwin>
         <gdrx_clutter_flag_filter>Off</gdrx_clutter_flag_filter>
         <cf_gip>On</cf_gip>
         <cf_gip_mode_const>On</cf_gip_mode_const>
         <cf_gip_max_width>0.5000</cf_gip_max_width>
         <cmd>On</cmd>
         <int_filter_phaseint>On</int_filter_phaseint>
         <int_filter_corr_mode>Interpolate</int_filter_corr_mode>
         <cf_dectree>Off</cf_dectree>
         <kdpautoscale>On</kdpautoscale>
         <triprecovery>Off</triprecovery>
         <midprf>0.00</midprf>
         <matchedfiltertype>SdfStaticLookupFilter</matchedfiltertype>
         <dynpw>0.95</dynpw>
         <dynpw_constantdc_enabled>On</dynpw_constantdc_enabled>
         <dynpw_constantdc>80</dynpw_constantdc>
         <dynpw_matchedrangeresolution>None</dynpw_matchedrangeresolution>
         <nyquistextensionmode>FreePRF</nyquistextensionmode>
         <nyquistextensionfactor>3.00</nyquistextensionfactor>
         <scanno>8</scanno>
         <scantotalno>12</scantotalno>
         <volno>28462</volno>
         <trackno>6139</trackno>
         <radarprefix_sp>262plc1_sq</radarprefix_sp>
         <radarprefix_rcp>262plc1_sq_lk</radarprefix_rcp>
         <acufaultstat>OK</acufaultstat>
         <gdrx5faultstat>OK</gdrx5faultstat>
         <gdrxtransmitfreq>5649.8447</gdrxtransmitfreq>
         <mpdbcufaultstat>OK</mpdbcufaultstat>
         <radargroundheight>93.3</radargroundheight>
         <radarLightningProtClass>0</radarLightningProtClass>
         <athattenuation>0.019</athattenuation>
         <rsplambda>5.31</rsplambda>
         <rsradarheight>122.3</rsradarheight>
         <rxfaultstat>OK</rxfaultstat>
         <spbantgain>45.38</spbantgain>
         <injectlossh>20.7</injectlossh>
         <spbdphtxcalpowkw>217.77</spbdphtxcalpowkw>
         <spbdpvantgain>45.37</spbdpvantgain>
         <injectlossv>21.08</injectlossv>
         <spbdpvrxloss>3.51</spbdpvrxloss>
         <spbdpvsysgainoffset>0.0</spbdpvsysgainoffset>
         <spbdpvtxcalpowkw>210.86</spbdpvtxcalpowkw>
         <spbdpvtxloss>0.66</spbdpvtxloss>
         <spbhorbeam>0.95</spbhorbeam>
         <spbphidpoffpw>0.00</spbphidpoffpw>
         <spbradomloss>0.2</spbradomloss>
         <spbrefindex>0.93</spbrefindex>
         <spbrxloss>3.32</spbrxloss>
         <spbsysgainoffset>0.0</spbsysgainoffset>
         <spbtxcalpowkw>217.77</spbtxcalpowkw>
         <spbtxloss>0.66</spbtxloss>
         <spbverbeam>0.99</spbverbeam>
         <txfaultstat>OK</txfaultstat>
         <txwarningstat>OK</txwarningstat>
         <rxbandwidth>1100000</rxbandwidth>
         <powerdiff>-0.15</powerdiff>
         <noise_power_dbz>-42.2651</noise_power_dbz>
         <noise_power_dbz_dpv>-42.3474</noise_power_dbz_dpv>
         <rspradconst>69.534</rspradconst>
         <rspdphradconst>69.534</rspdphradconst>
         <rspdpvradconst>69.884</rspdpvradconst>
         <spbtxcorr>4</spbtxcorr>
         <gdrxafcmode>20</gdrxafcmode>
         <gdrxafcfreq>0.0</gdrxafcfreq>
         <gdrxanctxfreq>59.8483</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:48:54.193" time="18:48:54" stopdatetimehighaccuracy="2024-10-05T18:49:14.192">
            <rayinfo rays="360" blobid="0" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="1" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="2" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="1">
         <posangle>1.5</posangle>
         <scanno>7</scanno>
         <trackno>6138</trackno>
         <gdrxanctxfreq>59.9227</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:48:33.191" time="18:48:33" stopdatetimehighaccuracy="2024-10-05T18:48:53.192">
            <rayinfo rays="360" blobid="3" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="4" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="5" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="2">
         <posangle>2.5</posangle>
         <scanno>6</scanno>
         <trackno>6137</trackno>
         <gdrxanctxfreq>60.0078</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:48:12.193" time="18:48:12" stopdatetimehighaccuracy="2024-10-05T18:48:32.192">
            <rayinfo rays="360" blobid="6" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="7" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="8" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="3">
         <posangle>3.5</posangle>
         <scanno>5</scanno>
         <trackno>6136</trackno>
         <noise_power_dbz_dpv>-42.3475</noise_power_dbz_dpv>
         <gdrxanctxfreq>60.0972</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:47:51.357" time="18:47:51" stopdatetimehighaccuracy="2024-10-05T18:48:11.355">
            <rayinfo rays="360" blobid="9" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="10" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="11" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="4">
         <dynv max="35.1257" min="-35.1257"/>
         <dynw max="11.6628" min="0.0457365"/>
         <posangle>5.4</posangle>
         <stoprange>170</stoprange>
         <highprf>882</highprf>
         <lowprf>661.50</lowprf>
         <dynpw>1.29</dynpw>
         <antspeed>16</antspeed>
         <int_filter_hconf>6</int_filter_hconf>
         <int_filter_vconf>6</int_filter_vconf>
         <timesamp>24</timesamp>
         <cf_gip_width>0.3</cf_gip_width>
         <scanno>4</scanno>
         <trackno>6135</trackno>
         <rxbandwidth>800000</rxbandwidth>
         <noise_power_dbz>-44.9218</noise_power_dbz>
         <noise_power_dbz_dpv>-45.0041</noise_power_dbz_dpv>
         <rspradconst>68.206</rspradconst>
         <rspdphradconst>68.206</rspdphradconst>
         <rspdpvradconst>68.556</rspdpvradconst>
         <gdrxanctxfreq>60.2026</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:47:27.754" time="18:47:27" stopdatetimehighaccuracy="2024-10-05T18:47:50.252">
            <rayinfo rays="360" blobid="12" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="13" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="340" type="dBZ" blobid="14" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="5">
         <posangle>9.1</posangle>
         <antspeed>20</antspeed>
         <timesamp>26</timesamp>
         <calibslottypes>DATA,NOISE</calibslottypes>
         <cf_gip_width>0.3</cf_gip_width>
         <scanno>3</scanno>
         <trackno>6134</trackno>
         <gdrxanctxfreq>60.2829</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:47:08.491" time="18:47:08" stopdatetimehighaccuracy="2024-10-05T18:47:26.494">
            <rayinfo rays="360" blobid="15" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="16" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="17" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="6">
         <posangle>15.0</posangle>
         <timesamp>24</timesamp>
         <calibslottypes>DATA,NOISE</calibslottypes>
         <antspeed>22</antspeed>
         <cf_gip_width>0.3</cf_gip_width>
         <int_filter_mode>2</int_filter_mode>
         <int_filter_hconf>16</int_filter_hconf>
         <int_filter_vconf>16</int_filter_vconf>
         <scanno>2</scanno>
         <trackno>6133</trackno>
         <gdrxanctxfreq>60.111</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:46:50.621" time="18:46:50" stopdatetimehighaccuracy="2024-10-05T18:47:06.984">
            <rayinfo rays="360" blobid="18" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="19" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="20" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
      <slice refid="7">
         <posangle>23.8</posangle>
         <cf_gip_width>0.3</cf_gip_width>
         <int_filter_mode>2</int_filter_mode>
         <int_filter_hconf>16</int_filter_hconf>
         <int_filter_vconf>16</int_filter_vconf>
         <antspeed>22</antspeed>
         <timesamp>24</timesamp>
         <calibslottypes>DATA,NOISE</calibslottypes>
         <scanno>1</scanno>
         <trackno>6132</trackno>
         <noise_power_dbz_dpv>-42.3475</noise_power_dbz_dpv>
         <gdrxanctxfreq>59.9191</gdrxanctxfreq>
         <slicedata date="2024-10-05" datetimehighaccuracy="2024-10-05T18:46:32.575" time="18:46:32" stopdatetimehighaccuracy="2024-10-05T18:46:48.938">
            <rayinfo rays="360" blobid="21" depth="16" refid="startangle"/>
            <rayinfo rays="360" blobid="22" depth="16" refid="stopangle"/>
            <rawdata rays="360" bins="250" type="dBZ" blobid="23" depth="8" max="95.5" min="-31.5"/>
         </slicedata>
      </slice>
   </scan>
   <sensorinfo id="LEG" type="gdrx" name="Legionowo">
      <lon>20.961110</lon>
      <lat>52.405250</lat>
      <alt>122.300000</alt>
      <wavelen>0.0531</wavelen>
      <beamwidth>0.95</beamwidth>
   </sensorinfo>
</volume>
<!-- END XML -->
<BLOB blobid="0" size="735" compression="qt">

So, basically, here I'm attemping to create a Python script that will perform next tasks:
1. Read RAINBOW volume dBZ radar data (path to the data file: dBZ.vol)
2. Reproject the radar data file from original raw polar coordinates that are in the file to cartesian coordinates, and later to EPSG:2180 coordinate system
3. Clip the radar data to the extent of Masovian region
4. Visualize and plot the processed radar data in a simple way
Though currently, unfortunately, I'm encountering with several issues in the Python code. Your task - is to fix all of them and provide a complete updated code. The issues are:
1. Color bar on the plot is being duplicated
2. The extent of radar data on the plot is way smaller than it should be (the original raw extent of the radar scanning is perfectly 125km from the radar site, a circle-like extent common for the weather radars + of course this should be accounted to the data cutting for specific masovian region extent), here's a basic function that might help with getting the masovian region extent:
def get_masovian_boundary():
    url = "https://raw.githubusercontent.com/ppatrzyk/polska-geojson/master/wojewodztwa/wojewodztwa-medium.geojson"
    # Fetch the content using requests
    response = requests.get(url)
    if response.status_code == 200:
        data = response.json()
        # Set the CRS to EPSG:4326 when creating the GeoDataFrame
        gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4326")
        masovian = gdf[gdf['nazwa'] == 'mazowieckie']
        return masovian
    else:
        print(f"Error fetching GeoJSON data: {response.status_code}")
        return None
But of course in our code we already have a data file - the masovian_boundary.geojson. But it's structure is literally the same 1:1.
3. The radar data plotted on the plot differs from my expectations. It looks like a small circle filled with dark blue color and random, grainy data not representing actual radar data storaged in the dBZ.vol that we need to decode and visualize

Here's also a console log returning when running the code that might significantly help you:
2024-10-05 21:36:01,864 - INFO - Parsing the .vol file...
2024-10-05 21:36:01,866 - INFO - Slice refid 0: rays 360, bins 250.
2024-10-05 21:36:01,866 - INFO - Slice refid 1: rays 360, bins 250.
2024-10-05 21:36:01,866 - INFO - Slice refid 2: rays 360, bins 250.
2024-10-05 21:36:01,867 - INFO - Slice refid 3: rays 360, bins 250.
2024-10-05 21:36:01,867 - INFO - Slice refid 4: rays 360, bins 340.
2024-10-05 21:36:01,867 - INFO - Slice refid 5: rays 360, bins 250.
2024-10-05 21:36:01,867 - INFO - Slice refid 6: rays 360, bins 250.
2024-10-05 21:36:01,868 - INFO - Slice refid 7: rays 360, bins 250.
2024-10-05 21:36:01,868 - INFO - Total slices with dBZ data: 8.

2024-10-05 21:36:01,868 - INFO - Extracting dBZ data from binary content...
2024-10-05 21:36:01,869 - WARNING - Slice 1 refid 0: Insufficient binary data. Expected 360000 bytes, got 182926 bytes. Adjusting number of bins.
2024-10-05 21:36:01,869 - INFO - Slice 1 refid 0: Extracted dBZ data with shape (360, 127).
2024-10-05 21:36:01,869 - WARNING - Slice 2 refid 1: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,869 - ERROR - Slice 2 refid 1: No data available. Skipping this slice.
2024-10-05 21:36:01,869 - WARNING - Slice 3 refid 2: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,869 - ERROR - Slice 3 refid 2: No data available. Skipping this slice.
2024-10-05 21:36:01,869 - WARNING - Slice 4 refid 3: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,870 - ERROR - Slice 4 refid 3: No data available. Skipping this slice.
2024-10-05 21:36:01,870 - WARNING - Slice 5 refid 4: Insufficient binary data. Expected 489600 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,870 - ERROR - Slice 5 refid 4: No data available. Skipping this slice.
2024-10-05 21:36:01,870 - WARNING - Slice 6 refid 5: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,871 - ERROR - Slice 6 refid 5: No data available. Skipping this slice.
2024-10-05 21:36:01,871 - WARNING - Slice 7 refid 6: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,871 - ERROR - Slice 7 refid 6: No data available. Skipping this slice.
2024-10-05 21:36:01,871 - WARNING - Slice 8 refid 7: Insufficient binary data. Expected 360000 bytes, got 46 bytes. Adjusting number of bins.
2024-10-05 21:36:01,871 - ERROR - Slice 8 refid 7: No data available. Skipping this slice.
2024-10-05 21:36:01,872 - INFO - Maximum number of bins across slices: 127
2024-10-05 21:36:01,872 - INFO - Combined dBZ data shape after padding and stacking: (360, 127)

2024-10-05 21:36:01,872 - INFO - Converting polar coordinates to Cartesian coordinates...
2024-10-05 21:36:01,874 - INFO - Coordinate conversion completed.

2024-10-05 21:36:01,874 - INFO - Creating GeoDataFrame...
2024-10-05 21:36:03,113 - INFO - GeoDataFrame creation completed.

2024-10-05 21:36:03,120 - INFO - Reprojecting GeoDataFrame to EPSG:2180...
2024-10-05 21:36:04,891 - INFO - Reprojection completed.

2024-10-05 21:36:05,452 - INFO - Loaded Masovian boundary successfully.

2024-10-05 21:36:05,560 - INFO - Reprojected Masovian boundary to EPSG:2180.

2024-10-05 21:36:05,576 - INFO - Clipping radar data to Masovian boundary...
2024-10-05 21:39:18,682 - INFO - Clipped data contains 45720 points.

2024-10-05 21:39:18,682 - INFO - Visualizing radar data...
2024-10-05 21:39:20,075 - INFO - Radar data visualization saved to output/processed_radar_data.png.

2024-10-05 21:39:20,076 - INFO - Radar data processing and visualization completed successfully.