Convert CHIRPS GeoTIFF in a folder to single NetCDF file with time dimension enabled that is CF-Compliant

GeoTIFF_to_NetCDF.py

#!/usr/bin/env python
"""
-------------------------------------------------------------------------------------------------------------
Convert raster GeoTIFF in a folder to single NetCDF file with time dimension enabled that is CF-Compliant
http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html
 
Based on Rich Signell's answer on StackExchange: https://gis.stackexchange.com/a/70487
This script was tested using TerraClimate monthly data. Adjustment is needed if using other timesteps data 
NCO (http://nco.sourceforge.net), GDAL, NetCDF4, numpy must be installed before using this script
 
Modified by
Benny Istanto, 
   UN World Food Programme, [email protected] [1]
   GOST/DECAT/DECDG, The World Bank, [email protected] [2]

Log:
[1] 2020-07-05 First draft
[2] 2023-07-15 Prevent the output shifted half-pixel
-------------------------------------------------------------------------------------------------------------
"""
import numpy as np
import datetime as dt
import os
from osgeo import gdal
gdal.UseExceptions()
import netCDF4
import re

ds = gdal.Open('../terraclimate/q/geotiff_idn/idn_cli_terraclimate_q_19580101.tif') # Data location

a = ds.ReadAsArray()
nlat,nlon = np.shape(a)

b = ds.GetGeoTransform() #bbox, interval
lon = np.arange(nlon) * b[1] + b[0] + (b[1]/2)  # add half the x pixel size to the lon
lat = np.arange(nlat) * b[5] + b[3] + (b[5]/2)  # add half the y pixel size to the lat
lat = np.flipud(lat)  # flip the latitudes

basedate = dt.datetime(1958,1,1,0,0,0)

# Create NetCDF file
nco = netCDF4.Dataset('../terraclimate/q/nc_idn/idn_cli_terraclimate_q_1958_2022.nc','w',clobber = True) # Output name

# Create dimensions, variables and attributes:
nco.createDimension('lon',nlon)
nco.createDimension('lat',nlat)
nco.createDimension('time',None)

timeo = nco.createVariable('time','f4',('time'))
timeo.units = 'days since 1958-01-01 00:00:00'
timeo.standard_name = 'time'
timeo.calendar = 'gregorian'
timeo.axis = 'T'

lono = nco.createVariable('lon','f4',('lon'))
lono.units = 'degrees_east'
lono.standard_name = 'longitude'
lono.long_name = 'longitude'
lono.axis = 'X'

lato = nco.createVariable('lat','f4',('lat'))
lato.units = 'degrees_north'
lato.standard_name = 'latitude'
lato.long_name = 'latitude'
lato.axis = 'Y'

# Create container variable for CRS: lon/lat WGS84 datum
crso = nco.createVariable('crs','i4')
crso.long_name = 'Lon/Lat Coords in WGS84'
crso.grid_mapping_name='latitude_longitude'
crso.longitude_of_prime_meridian = 0.0
crso.semi_major_axis = 6378137.0
crso.inverse_flattening = 298.257223563

# Create float variable for precipitation data, with chunking
pcpo = nco.createVariable('q', 'f4',  ('time', 'lat', 'lon'),zlib=True,fill_value=-9999.)
pcpo.units = 'mm'
pcpo.standard_name = 'runoff_amount'
pcpo.long_name = 'runoff_amount'
pcpo.time_step = 'monthly'
pcpo.missing_value = -9999.
pcpo.geospatial_lat_min = -11.0416667
pcpo.geospatial_lat_max = 5.8750000
pcpo.geospatial_lon_min = 95.0416667
pcpo.geospatial_lon_max = 141.0416667
pcpo.grid_mapping = 'crs'
pcpo.set_auto_maskandscale(False)

# Additional attributes
nco.Conventions = 'CF-1.6'
nco.cdm_data_type = "GRID"
nco.title = "TerraClimate: monthly climate and climatic water balance for global land surfaces"
nco.summary = "This archive contains a dataset of high-spatial resolution (1/24°, ~4-km) monthly climate and climatic water balance for global terrestrial surfaces from 1958-2015. These data were created by using climatically aided interpolation, combining high-spatial resolution climatological normals from the WorldClim version 1.4 and version 2 datasets, with coarser resolution time varying (i.e. monthly) data from CRU Ts4.0 and JRA-55 to produce a monthly dataset of precipitation, maximum and minimum temperature, wind speed, vapor pressure, and solar radiation. TerraClimate additionally produces monthly surface water balance datasets using a water balance model that incorporates reference evapotranspiration, precipitation, temperature, and interpolated plant extractable soil water capacity."
nco.references = "Abatzoglou, J.T., S.Z. Dobrowski, S.A. Parks, and K.C. Hegewisch, 2017, High-resolution global dataset of monthly climate and climatic water balance from 1958-2015, submitted to Scientific Data."
nco.version = "v1.0"
nco.comments = "time variable denotes the first day of the given m."
nco.website = "http://climate.nkn.uidaho.edu/TerraClimate"
nco.date_created = "2023-10-09"
nco.creator_name = "Benny Istanto"
nco.creator_role = "Climate Geographer"
nco.creator_email = "[email protected]"
nco.institution = "The World Bank"
nco.note = "The data is developed to support regular updating procedure for drought analysis. This activities will support the World Bank to assess extreme dry and wet periods"

# Write lon,lat
lono[:] = lon
lato[:] = lat

pat = re.compile('idn_cli_terraclimate_q_([0-9]{8})\.tif$')
itime = 0

# Step through data, writing time and data to NetCDF
for root, dirs, files in os.walk('../terraclimate/q/geotiff_idn/'):
    dirs.sort()
    files.sort()
    for f in files:
        match = re.match(pat, f)
        if match:
            date_str = match.group(1)  # This will give you '19810101'
            date = dt.datetime.strptime(date_str, '%Y%m%d')  # Convert to datetime
            print(date)
            dtime = (date - basedate).total_seconds()/86400.
            timeo[itime] = dtime
            # precipitation
            pcp_path = os.path.join(root,f)
            print(pcp_path)
            pcp = gdal.Open(pcp_path)
            a = pcp.ReadAsArray() # data
            pcpo[itime,:,:] = np.flipud(a)  # flip the data in y-direction
            itime = itime + 1

nco.close()

print('Completed!')