extracting WRF data

[kcmonaka]

Good Day, I am trying to validate rainfall observed vs WRF model simulated. now I am trying to extract data corresponding to the station data available. I am using the "Haversine Distance Matching" am i doing it correctly?

I have already processed precipitation into daily data and I am using script below to extract the value of RFE from the WRF output


import xarray as xr
import pandas as pd
import numpy as np

# --- Function to calculate haversine distance between two points ---
def haversine(lat1, lon1, lat2, lon2):
"""
Calculate the great-circle distance (in kilometers) between two points
given in decimal degrees (arrays supported).
"""
R = 6371.0 # Radius of Earth in kilometers
lat1_rad = np.radians(lat1)
lat2_rad = np.radians(lat2)
delta_lat = np.radians(lat2 - lat1)
delta_lon = np.radians(lon2 - lon1)

a = np.sin(delta_lat / 2.0) ** 2 + np.cos(lat1_rad) * np.cos(lat2_rad) * np.sin(delta_lon / 2.0) ** 2
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))

return R * c # Distance in kilometers

# --- Define the station information ---
stations = [
{"name": "Moss` s Farm (Molalatau)", "lat": -22.07346, "lon": 28.59792},
{"name": "Maunatlala councillors Office", "lat": -22.5967099016421, "lon": 27.6378049864268},
{"name": "Selibe Phikwe Met Station", "lat": -21.97, "lon": 27.83},
{"name": "Plaatjan Border", "lat": -22.48525, "lon": 28.831611111111112},
{"name": "Mogapinyana Primary School", "lat": -22.3756944444444, "lon": 27.592277777777777},
{"name": "Semolale Primary School", "lat": -21.8761666666667, "lon": 28.830222222222222},
{"name": "Maunatala Police Station", "lat": -22.591951484229, "lon": 27.6384203891164},
{"name": "Tebogo Primary School (Selibe Phikwe)", "lat": -21.968056937257, "lon": 27.8391465830506},
{"name": "Nxai Pan Wildlife Camp", "lat": -19.8555249507367, "lon": 24.7822121192866},
{"name": "Plateau Primary School", "lat": -17.8005555555556, "lon": 25.161583333333333},
{"name": "Chobe Junior Secondary School", "lat": -17.78564, "lon": 25.16666},
{"name": "P.G Matante Airport Met Office", "lat": -21.15, "lon": 27.483333333333334},
{"name": "Mmashoro Primary School", "lat": -21.866666666666667, "lon": 26.433333333333334},
{"name": "Letlhakane Met Station", "lat": -21.42, "lon": 25.62},
{"name": "Orapa Police Station", "lat": -21.3265555555556, "lon": 25.370527777777777},
{"name": "Selibe Phikwe Prison", "lat": -21.9076452249287, "lon": 27.8675377437635},
{"name": "Kumakwane Primary School", "lat": -24.6642793775725, "lon": 25.6998079513267},
{"name": "Kgomokasitwa Primary School", "lat": -25.0904087039343, "lon": 25.6297727956324},
{"name": "Seaseole Primary School (Letlhakane)", "lat": -21.39318884658, "lon": 25.5826218739518},
{"name": "Masunga Police Station", "lat": -20.633333333333333, "lon": 27.45194444},
{"name": "Letlhakane Prison", "lat": -21.4066285969262, "lon": 25.5898117598271},
{"name": "Tobane Crops", "lat": -21.9110285374283, "lon": 28.0416840644143},
{"name": "Maun Airport Met Office", "lat": -19.9764444444444, "lon": 23.42688888888889},
{"name": "Senete Primary School", "lat": -20.3224722222222, "lon": 27.131222222222224},
{"name": "Dagwi Primary School", "lat": -20.26175, "lon": 27.25},
{"name": "Mogoditshane Police Station", "lat": -24.6073764976606, "lon": 25.8550043089964},
{"name": "Pandamatenga Met Station", "lat": -18.5445, "lon": 25.63577777777778},
{"name": "Motswakae`s Farm (Tonota)", "lat": -21.43783, "lon": 27.46338},
{"name": "Ditladi Primary School", "lat": -21.4618611111111, "lon": 27.53322222222222},
{"name": "Kasane Airport Met Office", "lat": -17.8294722222222, "lon": 25.16486111111111},
{"name": "Sechele Primary School", "lat": -20.7541944444444, "lon": 27.357277777777778},
{"name": "L.O Dialwa`s Farm (Lerala)", "lat": -22.78897, "lon": 27.76186},
{"name": "Mogoditshane Kgotla ", "lat": -24.626840316893, "lon": 25.8720607306874},
{"name": "Masunga Primary School", "lat": -20.6203055555556, "lon": 27.445972222222224},
{"name": "Tshesebe Police Station", "lat": -20.7523888888889, "lon": 27.600666666666665},
{"name": "Kazungula Primary School", "lat": -17.7955, "lon": 25.234472222222223},
{"name": "Kazungula Police Station", "lat": -17.8041388888889, "lon": 25.24522222222222},
{"name": "Thamaga Police", "lat": -24.666666666666668, "lon": 25.533333333333335},
{"name": "Tirelo Primary School (Dibete)", "lat": -23.74812, "lon": 26.47174},
{"name": "St Joseph`s Primary School (Kgale)", "lat": -24.714190145595, "lon": 25.8717856693041},
{"name": "Moiyabana Primary School", "lat": -22.38, "lon": 26.24},
{"name": "Mabalane Primary School", "lat": -24.6226, "lon": 26.37483},
{"name": "Mandunyane Primary School", "lat": -21.3841388888889, "lon": 27.415305555555555},
{"name": "Kazungula Road Border", "lat": -17.8151944444444, "lon": 25.261666666666667},
{"name": "Ithuteng Junior Secondary School (Mochudi)", "lat": -24.3499166666667, "lon": 26.151027777777777},
{"name": "Sir Seretse Khama Airport Met Office", "lat": -24.7, "lon": 25.916666666666668},
{"name": "Motshegaletau Primary School", "lat": -22.34, "lon": 26.33},
{"name": "Oliphants Drift Police Station", "lat": -24.1925, "lon": 26.86513888888889},
{"name": "Artesia Primary School", "lat": -24.01, "lon": 26.317166666666665},
{"name": "Gaborone Village Met HQ", "lat": -24.666666666666668, "lon": 25.916666666666668},
{"name": "Sowa Met Station", "lat": -20.633333333333333, "lon": 25.983333333333334},
{"name": "Palapye Police Station", "lat": -22.509230701939, "lon": 27.1103260965648},
{"name": "Shakawe Met Station", "lat": -18.366666666666667, "lon": 21.833333333333332},
{"name": "Swaneng Hill School (Serowe)", "lat": -22.4140174599367, "lon": 26.7473119782236},
{"name": "Mosetse Primary School", "lat": -20.6495277777778, "lon": 26.64988888888889},
{"name": "Gammathari`s Farm (Kgope)", "lat": -24.294825, "lon": 25.9641683},
{"name": "Mogome Primary School", "lat": -22.7827438239323, "lon": 26.8977230844172},
{"name": "Borwa Junior Secondary School (Bokaa)", "lat": -24.434306, "lon": 26.044528},
{"name": "Mookane Primary School", "lat": -23.6841118764018, "lon": 26.659938832241},
{"name": "Morama Junior Secondary School (Jwaneng)", "lat": -24.60319, "lon": 24.73664},
{"name": "Jwaneng Met Station", "lat": -24.58, "lon": 24.8},
{"name": "Mahalapye Met Station", "lat": -23.116666666666667, "lon": 26.833333333333332},

]

# --- Load the WRF NetCDF file ---
nc_file = "wrfout_d02_2023-12-26_daily_precip.nc"
ds = xr.open_dataset(nc_file)

# --- Check required variables ---
if "XLAT" not in ds.variables or "XLONG" not in ds.variables:
raise KeyError("Variables 'XLAT' or 'XLONG' are missing in the NetCDF file.")

# Extract 2D latitude and longitude arrays
lats = ds["XLAT"].values.squeeze()
lons = ds["XLONG"].values.squeeze()

# --- Load precipitation ---
precip_var = "TOTAL_PRECIP" # Adjust if your precipitation variable name is different
if precip_var not in ds.variables:
raise KeyError(f"Variable '{precip_var}' not found in the NetCDF file.")

precip = ds[precip_var].values.squeeze()

# --- Prepare list for storing extracted station data ---
station_data = []

# --- Process each station ---
for station in stations:
name = station["name"]
lat = station["lat"]
lon = station["lon"]

# Calculate haversine distance between the station and all grid points
distance = haversine(lat, lon, lats, lons)

# Find the grid point with the minimum distance
nearest_lat_idx, nearest_lon_idx = np.unravel_index(distance.argmin(), distance.shape)

# Extract precipitation at the nearest grid point
precip_value = precip[nearest_lat_idx, nearest_lon_idx]

# Append station data
station_data.append({
"Station": name,
"Latitude": lat,
"Longitude": lon,
"Precipitation (mm)": precip_value
})

# --- Convert to DataFrame and save ---
df = pd.DataFrame(station_data)

# Save to CSV
output_csv = "station_precip_haversine.csv"
df.to_csv(output_csv, index=False)

print(f" Station data saved to {output_csv}")

 

[Ming Chen]

Hi,

I take a quick look at your script and it looks fine to me. Have you checked the results and whether they look reasonable?

Just in case, can you also try the function ll_to_xy provided in wrf-python package? Please see details here:

wrf.ll_to_xy — wrf-python 1.3.4.1 documentation

I won't expect identical results using the two different methods, However, I suppose the results should be at least close. Please let me know if you have any update. Thanks in advance.

 

[William.Hatheway]

I'm going to include a few snippets of my python code that might help you.

# Check if the correct arguments were provided
if len(sys.argv) not in [4, 6]:
    print(
        "\nEnter the required arguments: path_wrf, domain, city, [latitude and longitude (optional)]\n"
        "For example: script_name.py /home/WRF/test/em_real d01 Paris [48.8566 2.3522]\n"
    )
    sys.exit()

# Define the path where the netcdf files are and the domain to be used
path_wrf = sys.argv[1]
domain = sys.argv[2]
city = sys.argv[3]

# Define latitude and longitude
if len(sys.argv) == 6:
    try:
        lat = float(sys.argv[4])
        lon = float(sys.argv[5])
    except ValueError:
        print(
            "Invalid latitude or longitude. Please enter valid decimal numbers for latitude and longitude."
        )
        sys.exit()
else:
    # Prompt user for GPS coordinates if not provided via command-line
    while True:
        try:
            lat = float(
                input("Enter latitude in decimal format (e.g., 48.8566): ").strip()
            )
            break
        except ValueError:
            print(
                "Invalid input. Please enter latitude in decimal format (e.g., 48.8566)."
            )

    while True:
        try:
            lon = float(
                input("Enter longitude in decimal format (e.g., 2.3522): ").strip()
            )
            break
        except ValueError:
            print(
                "Invalid input. Please enter longitude in decimal format (e.g., 2.3522)."
            )


    # Get latitude and longtiude locations for Skew-T
    lat_lon = [lat, lon]
    x_y = wrf.ll_to_xy(ncfile, lat_lon[0], lat_lon[1])


    T = T1[:, x_y[1], x_y[0]] * units.degC

hopefully this can help you