Georeferencing downloaded CMIP6 data and extracting an area of interest

CMIP6 data is publicly available and accessible through various data repositories and portals, such as the Earth System Grid Federation (ESGF). This article completes the first one I wrote downloading CMIP6 data. You can access it here.

CMIP6 data are not typically georeferenced in the traditional sense for direct use in GIS application. CMIP6 data use a 0 to 360-degree longitude range because it aligns with the 360-degree system used in many global climate models. This system ensures that the longitudinal data wraps around the globe seamlessly, with 0 degrees (or 360 degrees) corresponding to the Prime Meridian (Greenwich Meridian) and 180 degrees representing the International Date Line. This convention is particularly useful for global climate simulations and model intercomparisons, as it simplifies the handling of data at the boundaries of the longitude range. It also helps in maintaining consistency across different models and datasets.

This georeferencing system, can be a barrier to work with CMIP6 data especially in GIS application. It is therefore, necessary to do some preprocessing. In this article, I will show you two important preprocessing steps:

• Georeferencing CMIP6 data: By georeferencing here I mean changing longitude from 0–360 degree system to -180 to 180 degree
• Extract an area of interest (AOI) by providing a bounding box.

The script associated with the article can be downloaded here.

Without further due, below are the necessary steps.

1. Importing relevant libraries

import xarray as xr
import numpy as np

2. Creating a function for georeferencing and extracting the AOI

def georef_clip (input_file, aoi, output_file): 
    dataset = xr.open_dataset(input_file) # Opening the CMIP6 dataset

    # Georeferencing step
    dataset["lon"] = np.where(dataset["lon"] > 180, dataset["lon"] - 360, dataset["lon"]) # Changing the 0-360 longitude format to use the -180 to 180 format
    dataset = dataset.sortby(dataset.lon) # Sorting the data by longitude values

    # Extracting the area of interest (AOI)

        # Finding the coordinates in the dataset that is closest to the provided AOI. 
    lat_max = np.argmin(np.abs(dataset["lat"].values - aoi[0]))
    lat_min = np.argmin(np.abs(dataset["lat"].values - aoi[1]))
    lon_max = np.argmin(np.abs(dataset["lon"].values - aoi[2]))
    lon_min = np.argmin(np.abs(dataset["lon"].values - aoi[3]))

        # Getting the id of the closest coordinates
    lon_min_id = dataset["lon"].values[lon_min]
    lon_max_id = dataset["lon"].values[lon_max]
    lat_min_id = dataset["lat"].values[lat_min]
    lat_max_id = dataset["lat"].values[lat_max]

        # Clipping the data using the identified coordinates in the dataset that close to the provided AOI
    clipped_data = dataset.sel(lon=slice(lon_min_id, lon_max_id), lat=slice(lat_min_id, lat_max_id))

    # Saving the result
    clipped_data.to_netcdf(output_file)

    print (f"Process completed. Data saved here {output_file}")

3. Launching the georeferencing process and extracting the AOI

# Setting the parameters to be used by the function
input_file = "C:/Users/ilung/Documents/uas_day_IPSL-CM6A-LR_ssp245_r1i1p1f1_gr_20150101-21001231.nc" # input CMIP6 file
aoi = [15, -3, 3, -14] # The AOI bounding box. Must use the following format [lat_max, lat_min, lon_max, lon_min]
output_file = "C:/Users/ilung/Documents/uas_day_clip.nc" # Output file

# Launching the function to georeference and extract the AOI
georef_clip(input_file, aoi, output_file)