Last active
November 11, 2025 15:07
-
-
Save planemad/00e99b8bd6e6464dd6ee0fc6aa79aacf to your computer and use it in GitHub Desktop.
India Meteorological Department (IMD) high resolution gridded rainfall data convertor for GIS
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| import os | |
| def grd_to_asc(input_grd, output_tif, year=None): | |
| """ | |
| Convert high resolution 0.25 x 0.25 deg IMD daily .grd file (IMD4) to Arc/Info ASCII Grid .asc for QGIS | |
| Download data from https://imdpune.gov.in/lrfindex.php | |
| EXAMPLES | |
| # Get annual total | |
| python3 imd4_grd_to_asc.py Rainfall_ind2024_rfp25.grd | |
| # Get monthly totals (12 files: Jan, Feb, Mar, etc.) | |
| python3 imd4_grd_to_asc.py Rainfall_ind2024_rfp25.grd --monthly | |
| # Get all 366 daily files | |
| python3 imd4_grd_to_asc.py Rainfall_ind2024_rfp25.grd --daily | |
| Based on IMD4 specification: | |
| - Data arranged as 135x129 grid points (longitude x latitude) | |
| - First data: 6.5°N & 66.5°E | |
| - Last data: 38.5°N & 100.0°E | |
| - For each day: data stored as ((RF(day,lon,lat),lon=1,135),lat=1,129) | |
| - 365/366 records for non-leap/leap years | |
| References: | |
| - IMD4 Paper https://mausamjournal.imd.gov.in/index.php/MAUSAM/article/view/851 | |
| - https://github.com/answerquest/IMD-grid-data-work | |
| - https://github.com/Manisht9/imdData/blob/main/imd_data.py | |
| """ | |
| # IMD grid specifications | |
| nlon = 135 # Longitude points: 66.5°E to 100.0°E | |
| nlat = 129 # Latitude points: 6.5°N to 38.5°N | |
| lon_start = 66.5 | |
| lat_start = 6.5 | |
| resolution = 0.25 | |
| # Calculate bounds | |
| lon_end = lon_start + (nlon - 1) * resolution # 100.0°E | |
| lat_end = lat_start + (nlat - 1) * resolution # 38.5°N | |
| # Determine if leap year | |
| if year is None: | |
| import re | |
| match = re.search(r'(\d{4})', os.path.basename(input_grd)) | |
| if match: | |
| year = int(match.group(1)) | |
| is_leap = False | |
| if year: | |
| is_leap = (year % 4 == 0 and year % 100 != 0) or (year % 400 == 0) | |
| ndays = 366 if is_leap else 365 | |
| print(f"Year: {year if year else 'Unknown'}") | |
| print(f"Leap year: {'Yes' if is_leap else 'No'}") | |
| print(f"Days: {ndays}") | |
| print(f"Grid: {nlon} lon × {nlat} lat") | |
| print(f"Lon: {lon_start}°E to {lon_end}°E") | |
| print(f"Lat: {lat_start}°N to {lat_end}°N") | |
| # Read binary data (little-endian 32-bit float) | |
| with open(input_grd, 'rb') as f: | |
| data = np.fromfile(f, dtype='<f4') | |
| total_points = nlon * nlat | |
| expected_size = ndays * total_points | |
| print(f"\nData points read: {len(data):,}") | |
| print(f"Expected: {ndays} × {nlon} × {nlat} = {expected_size:,}") | |
| if len(data) != expected_size: | |
| print(f"⚠ Warning: Size mismatch!") | |
| actual_days = len(data) // total_points | |
| print(f"Using {actual_days} days based on file size") | |
| ndays = actual_days | |
| # Reshape according to IMD format: (days, lon, lat) | |
| # Each day contains: for each lat (S to N), all lon values (W to E) | |
| # So data is actually stored as (days, lat, lon) in the file | |
| data = data.reshape(ndays, nlat, nlon) | |
| # Now data[day, lat, lon] where: | |
| # - lat=0 corresponds to 6.5°N (south) | |
| # - lat=128 corresponds to 38.5°N (north) | |
| # - lon=0 corresponds to 66.5°E (west) | |
| # - lon=134 corresponds to 100.0°E (east) | |
| print(f"\nCalculating annual total...") | |
| # Sum across all days to get annual total | |
| annual_total = np.sum(data, axis=0) # Result: (nlat, nlon) | |
| # Replace missing/invalid values with NaN | |
| annual_total[annual_total < -100] = np.nan | |
| annual_total[annual_total > 20000] = np.nan # Sanity check | |
| # Get statistics | |
| valid_data = annual_total[~np.isnan(annual_total)] | |
| print(f"\nAnnual rainfall statistics:") | |
| print(f" Min: {np.nanmin(annual_total):.2f} mm") | |
| print(f" Max: {np.nanmax(annual_total):.2f} mm") | |
| print(f" Mean: {np.nanmean(annual_total):.2f} mm") | |
| print(f" Valid points: {len(valid_data):,} / {annual_total.size:,}") | |
| # Create ASCII Grid | |
| create_ascii_grid(annual_total, output_tif.replace('.tif', '.asc'), | |
| lon_start, lat_start, resolution, nlat, nlon) | |
| return annual_total, data | |
| def create_ascii_grid(data, output_asc, xll, yll, cellsize, nrows, ncols): | |
| """ | |
| Create ESRI ASCII Grid format. | |
| Data should be (nlat, nlon) where: | |
| - data[0,:] is the southernmost latitude | |
| - data[-1,:] is the northernmost latitude | |
| ASCII Grid format expects data from NORTH to SOUTH. | |
| """ | |
| with open(output_asc, 'w') as f: | |
| # Write header | |
| f.write(f"ncols {ncols}\n") | |
| f.write(f"nrows {nrows}\n") | |
| f.write(f"xllcorner {xll}\n") | |
| f.write(f"yllcorner {yll}\n") | |
| f.write(f"cellsize {cellsize}\n") | |
| f.write(f"NODATA_value -9999\n") | |
| # Write data from north to south | |
| # data[0] is south (6.5N), data[-1] is north (38.5N) | |
| # So we write in reverse order: data[-1], data[-2], ..., data[0] | |
| for lat_idx in range(nrows-1, -1, -1): | |
| row_data = data[lat_idx, :] | |
| row_str = ' '.join([f"{val:.2f}" if not np.isnan(val) else "-9999" | |
| for val in row_data]) | |
| f.write(row_str + '\n') | |
| print(f"\n✓ ASCII grid saved: {output_asc}") | |
| def export_daily_layers(data, base_name, lon_start, lat_start, resolution, nlat, nlon): | |
| """ | |
| Export each day as a separate ASCII grid file | |
| """ | |
| ndays = data.shape[0] | |
| output_dir = f"{base_name}_daily" | |
| os.makedirs(output_dir, exist_ok=True) | |
| print(f"\nExporting {ndays} daily layers to {output_dir}/") | |
| for day in range(ndays): | |
| day_data = data[day, :, :] # (nlat, nlon) | |
| day_data[day_data < -100] = np.nan | |
| day_data[day_data > 1000] = np.nan # Sanity check for daily rainfall | |
| output_file = os.path.join(output_dir, f"day_{day+1:03d}.asc") | |
| create_ascii_grid(day_data, output_file, lon_start, lat_start, | |
| resolution, nlat, nlon) | |
| if (day + 1) % 50 == 0: | |
| print(f" Processed {day + 1}/{ndays} days...") | |
| print(f"✓ All {ndays} daily layers exported") | |
| def export_monthly_totals(data, base_name, lon_start, lat_start, resolution, nlat, nlon, year): | |
| """ | |
| Export monthly total rainfall | |
| """ | |
| is_leap = (year % 4 == 0 and year % 100 != 0) or (year % 400 == 0) | |
| # Days in each month | |
| if is_leap: | |
| days_in_month = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] | |
| else: | |
| days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] | |
| month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', | |
| 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] | |
| output_dir = f"{base_name}_monthly" | |
| os.makedirs(output_dir, exist_ok=True) | |
| print(f"\nExporting monthly totals to {output_dir}/") | |
| day_idx = 0 | |
| for month_idx, (month_name, ndays) in enumerate(zip(month_names, days_in_month)): | |
| # Sum days for this month | |
| month_data = np.sum(data[day_idx:day_idx+ndays, :, :], axis=0) | |
| month_data[month_data < -100] = np.nan | |
| output_file = os.path.join(output_dir, f"{month_idx+1:02d}_{month_name}.asc") | |
| create_ascii_grid(month_data, output_file, lon_start, lat_start, | |
| resolution, nlat, nlon) | |
| print(f" {month_name}: {np.nanmean(month_data):.1f} mm (avg)") | |
| day_idx += ndays | |
| print(f"✓ All 12 monthly layers exported") | |
| if __name__ == "__main__": | |
| import sys | |
| if len(sys.argv) < 2: | |
| print("Usage: python script.py input.grd [output.asc] [options]") | |
| print("\nExample: python script.py Rainfall_ind2024_rfp25.grd") | |
| print("\nOptions:") | |
| print(" --daily Export all 365/366 daily layers") | |
| print(" --monthly Export 12 monthly total layers") | |
| print("\nConverts IMD 0.25° gridded rainfall to ASCII Grid (QGIS compatible)") | |
| sys.exit(1) | |
| input_file = sys.argv[1] | |
| export_daily = '--daily' in sys.argv | |
| export_monthly = '--monthly' in sys.argv | |
| output_file = None | |
| for arg in sys.argv[2:]: | |
| if not arg.startswith('--'): | |
| output_file = arg | |
| break | |
| if output_file is None: | |
| output_file = input_file.replace('.grd', '_annual.asc') | |
| print(f"\n{'='*70}") | |
| print(f"IMD RAINFALL DATA CONVERTER") | |
| print(f"{'='*70}") | |
| print(f"Input: {input_file}") | |
| print(f"Output: {output_file}") | |
| print(f"{'='*70}\n") | |
| annual_total, daily_data = grd_to_asc(input_file, output_file) | |
| # Extract year for monthly export | |
| import re | |
| match = re.search(r'(\d{4})', os.path.basename(input_file)) | |
| year = int(match.group(1)) if match else 2024 | |
| if export_monthly: | |
| base_name = output_file.replace('.asc', '') | |
| export_monthly_totals(daily_data, base_name, 66.5, 6.5, 0.25, 129, 135, year) | |
| if export_daily: | |
| base_name = output_file.replace('.asc', '') | |
| export_daily_layers(daily_data, base_name, 66.5, 6.5, 0.25, 129, 135) | |
| print(f"\n{'='*70}") | |
| print("✓ CONVERSION COMPLETE") | |
| print(f"{'='*70}") | |
| print("\nTo open in QGIS:") | |
| print("1. Layer → Add Layer → Add Raster Layer") | |
| print(f"2. Select: {output_file}") | |
| print("3. The file contains ANNUAL TOTAL rainfall") | |
| print("4. CRS: EPSG:4326 (WGS 84)") | |
| print("\nExpected rainfall patterns:") | |
| print(" • High: Western Ghats (west coast), Northeast India") | |
| print(" • Medium: Eastern coast, Indo-Gangetic plains") | |
| print(" • Low: Northwest (Rajasthan), rain shadow areas") | |
| if export_monthly: | |
| print(f"\n Monthly totals: {output_file.replace('.asc', '')}_monthly/") | |
| if export_daily: | |
| print(f" Daily data: {output_file.replace('.asc', '')}_daily/") | |
| print(f"{'='*70}\n") |
Author
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Code credits: Claude Sonnet 4.5