|
| 1 | +# LIBRARIES ___________________________________________________________________ |
| 2 | +import os # For handling file paths |
| 3 | +import pandas as pd # For handling tabular data (VIC input/output) |
| 4 | +import xarray as xr # For working with NetCDF climate data |
| 5 | + |
| 6 | +# PRE-PROCESSING _____________________________________________________________ |
| 7 | + |
| 8 | +# Define the spatial domain you want to process (latitude and longitude range) |
| 9 | +min_lat, max_lat = 15, 35 |
| 10 | +min_lon, max_lon = 90, 105 |
| 11 | + |
| 12 | +# Define the range of years to process |
| 13 | +sta_year = 2022 |
| 14 | +end_year = 2023 |
| 15 | + |
| 16 | +# Set directories for input data and output VIC files |
| 17 | +main_dir = os.getcwd() |
| 18 | +prec_dir = os.path.join(main_dir, "prec") # Precipitation data folder |
| 19 | +tmax_dir = os.path.join(main_dir, "tmax") # Maximum temperature folder |
| 20 | +tmin_dir = os.path.join(main_dir, "tmin") # Minimum temperature folder |
| 21 | +wind_dir = os.path.join(main_dir, "wind") # Wind speed folder |
| 22 | +outp_dir = os.path.join(main_dir, "vic_inp") # Output folder for VIC files |
| 23 | + |
| 24 | +# Ensure the output directory exists |
| 25 | +os.makedirs(outp_dir, exist_ok=True) |
| 26 | + |
| 27 | +# PROCESSING __________________________________________________________________ |
| 28 | + |
| 29 | +# -------------------- Precipitation -------------------- |
| 30 | +print("Processing precipitation data: START") |
| 31 | +prec = [] # List to store each year's data |
| 32 | +years = range(sta_year, end_year+1) |
| 33 | +for year in years: |
| 34 | + print("year: "+str(year)) |
| 35 | + file_name = f"chirps-v2.0.{year}.days_p05.nc" # File naming convention |
| 36 | + file_path = os.path.join(prec_dir, file_name) |
| 37 | + |
| 38 | + # Open NetCDF file with xarray |
| 39 | + data = xr.open_dataset(file_path) |
| 40 | + |
| 41 | + # Slice the data to the defined spatial domain |
| 42 | + sliced_data = data.sel({"latitude": slice(min_lat, max_lat), |
| 43 | + "longitude": slice(min_lon, max_lon)}) |
| 44 | + prec.append(sliced_data) |
| 45 | + |
| 46 | +# Concatenate all years along the time dimension |
| 47 | +prec = xr.concat(prec, dim="time") |
| 48 | +print("Processing precipitation data: DONE") |
| 49 | + |
| 50 | +# -------------------- Maximum Temperature -------------------- |
| 51 | +print("Processing maximum temperature data: START") |
| 52 | +file_name = "era5.02m.tmax.1200.22-23.nc" |
| 53 | +file_path = os.path.join(tmax_dir, file_name) |
| 54 | + |
| 55 | +# Open NetCDF file |
| 56 | +data = xr.open_dataset(file_path) |
| 57 | + |
| 58 | +# ERA5 latitude may be from north to south, invert if needed |
| 59 | +inverted_data = data.isel({"latitude": slice(None, None, -1)}) |
| 60 | + |
| 61 | +# Slice to defined spatial domain |
| 62 | +tmax = inverted_data.sel({"latitude": slice(min_lat, max_lat), |
| 63 | + "longitude": slice(min_lon, max_lon)}) |
| 64 | +print("Processing maximum temperature data: DONE") |
| 65 | + |
| 66 | +# -------------------- Minimum Temperature -------------------- |
| 67 | +print("Processing minimum temperature data: START") |
| 68 | +file_name = "era5.02m.tmin.0000.22-23.nc" |
| 69 | +file_path = os.path.join(tmin_dir, file_name) |
| 70 | + |
| 71 | +data = xr.open_dataset(file_path) |
| 72 | +inverted_data = data.isel({"latitude": slice(None, None, -1)}) |
| 73 | +tmin = inverted_data.sel({"latitude": slice(min_lat, max_lat), |
| 74 | + "longitude": slice(min_lon, max_lon)}) |
| 75 | +print("Processing minimum temperature data: DONE") |
| 76 | + |
| 77 | +# -------------------- Wind Speed -------------------- |
| 78 | +print("Processing windspeed data: START") |
| 79 | +file_name = "era5.10m.wind.0600.22-23.nc" |
| 80 | +file_path = os.path.join(wind_dir, file_name) |
| 81 | + |
| 82 | +data = xr.open_dataset(file_path) |
| 83 | +inverted_data = data.isel({"latitude": slice(None, None, -1)}) |
| 84 | +wind = inverted_data.sel({"latitude": slice(min_lat, max_lat), |
| 85 | + "longitude": slice(min_lon, max_lon)}) |
| 86 | +print("Processing windspeed data: DONE") |
| 87 | + |
| 88 | +# EXTRACTING AND WRITING OUT DATA FOR VIC _____________________________________ |
| 89 | +print("Extracting and writing out data for VIC: START") |
| 90 | + |
| 91 | +# Load locations (longitude, latitude) where VIC input files are needed |
| 92 | +loc = pd.read_csv("loc_sample.txt", sep="\t", header=None, names=["lon", "lat"]) |
| 93 | + |
| 94 | +# Loop over each location |
| 95 | +for index, row in loc.iterrows(): |
| 96 | + lon = row['lon'] |
| 97 | + lat = row['lat'] |
| 98 | + print("lat: "+str(lat)+", lon: "+str(lon)) |
| 99 | + |
| 100 | + # Extract the nearest grid point data for each variable |
| 101 | + sel_prec = prec["precip"].sel(latitude=lat, longitude=lon, method="nearest") |
| 102 | + sel_tmax = tmax["mx2t"].sel(latitude=lat, longitude=lon, method="nearest") |
| 103 | + sel_tmin = tmin["mn2t"].sel(latitude=lat, longitude=lon, method="nearest") |
| 104 | + sel_wind = wind["fg10"].sel(latitude=lat, longitude=lon, method="nearest") |
| 105 | + |
| 106 | + # Create a pandas DataFrame with VIC-required variables |
| 107 | + df = pd.DataFrame({ |
| 108 | + "prec": sel_prec.values, # Precipitation (mm/day) |
| 109 | + "tmax": sel_tmax.values - 273.15, # Max temperature converted from K to °C |
| 110 | + "tmin": sel_tmin.values - 273.15, # Min temperature converted from K to °C |
| 111 | + "wind": sel_wind.values # Wind speed (m/s) |
| 112 | + }) |
| 113 | + |
| 114 | + # Define output filename and write to file (tab-separated, no header/index) |
| 115 | + filename = f"gf_{lon:.4f}_{lat:.4f}" # Original format, no .txt extension |
| 116 | + file_path = os.path.join(outp_dir, filename) |
| 117 | + df.to_csv(file_path, sep="\t", index=False, header=False, float_format="%.6f") |
| 118 | + |
| 119 | +print("Extracting and writing out data for VIC: DONE") |
0 commit comments