Convert WRF NetCDF to TIFF

Question

I have a NetCDF file, output of a WRF model, which contains a time series of some variables. You can download it here (Warning: 7GB download). Here is the gdalinfo output:

Driver: netCDF/Network Common Data Format
Files: wrfout_d03_Florence_Control.nc
       wrfout_d03_Florence_Control.nc.aux.xml
Size is 512, 512
Metadata:
  NC_GLOBAL#AERCU_FCT=1
  NC_GLOBAL#AERCU_OPT=0
  NC_GLOBAL#AER_ANGEXP_OPT=1
  NC_GLOBAL#AER_ANGEXP_VAL=1.3
  NC_GLOBAL#AER_AOD550_OPT=1
  NC_GLOBAL#AER_AOD550_VAL=0.12
  NC_GLOBAL#AER_ASY_OPT=1
  NC_GLOBAL#AER_ASY_VAL=0.89999998
  NC_GLOBAL#AER_OPT=0
  NC_GLOBAL#AER_SSA_OPT=1
  NC_GLOBAL#AER_SSA_VAL=0.85000002
  NC_GLOBAL#AER_TYPE=1
  NC_GLOBAL#AUTO_LEVELS_OPT=2
  NC_GLOBAL#BLDT=0
  NC_GLOBAL#BL_PBL_PHYSICS=2
  NC_GLOBAL#BOTTOM-TOP_GRID_DIMENSION=36
  NC_GLOBAL#BOTTOM-TOP_PATCH_END_STAG=36
  NC_GLOBAL#BOTTOM-TOP_PATCH_END_UNSTAG=35
  NC_GLOBAL#BOTTOM-TOP_PATCH_START_STAG=1
  NC_GLOBAL#BOTTOM-TOP_PATCH_START_UNSTAG=1
  NC_GLOBAL#BUCKET_J=1e+09
  NC_GLOBAL#BUCKET_MM=100
  NC_GLOBAL#CDI=Climate Data Interface version 2.0.4 (https://mpimet.mpg.de/cdi)
  NC_GLOBAL#CDO=Climate Data Operators version 2.0.4 (https://mpimet.mpg.de/cdo)
  NC_GLOBAL#CEN_LAT=43.890499
  NC_GLOBAL#CEN_LON=11.279205
  NC_GLOBAL#Conventions=CF-1.6
  NC_GLOBAL#CUDT=0
  NC_GLOBAL#CU_PHYSICS=0
  NC_GLOBAL#DAMPCOEF=0.15000001
  NC_GLOBAL#DAMP_OPT=3
  NC_GLOBAL#DFI_OPT=0
  NC_GLOBAL#DIFF_6TH_FACTOR=0.050000001
  NC_GLOBAL#DIFF_6TH_OPT=0
  NC_GLOBAL#DIFF_6TH_SLOPEOPT=0
  NC_GLOBAL#DIFF_6TH_THRESH=0.1
  NC_GLOBAL#DIFF_OPT=1
  NC_GLOBAL#DT=5
  NC_GLOBAL#DVEG=5
  NC_GLOBAL#DX=1000
  NC_GLOBAL#DY=1000
  NC_GLOBAL#DZBOT=50
  NC_GLOBAL#DZSTRETCH_S=1.3
  NC_GLOBAL#DZSTRETCH_U=1.1
  NC_GLOBAL#ETAC=0.2
  NC_GLOBAL#FEEDBACK=0
  NC_GLOBAL#GFDDA_END_H=0
  NC_GLOBAL#GFDDA_INTERVAL_M=0
  NC_GLOBAL#GMT=0
  NC_GLOBAL#GRAV_SETTLING=0
  NC_GLOBAL#GRIDTYPE=C
  NC_GLOBAL#GRID_FDDA=0
  NC_GLOBAL#GRID_ID=1
  NC_GLOBAL#GRID_SFDDA=0
  NC_GLOBAL#GWD_OPT=0
  NC_GLOBAL#history=Fri Mar 03 14:15:25 2023: cdo mergetime TEMPORARY/CTL_wrfout_d01_2018-12-29_00:00:00.nc TEMPORARY/CTL_wrfout_d01_2019-05-13_00:00:00.nc TEMPORARY/CTL_wrfout_d01_2019-07-12_00:00:00.nc TEMPORARY/CTL_wrfout_d01_2019-10-15_00:00:00.nc wrfout_d03_Florence_Control.nc
Wed Nov  2 09:06:53 2022: ncks -d bottom_top,0,0 -d Time,0,3317 -v o3,PM2_5_DRY,PM10,no2,T2,MEBIO_ISOP,SWDOWN,dvel_o3,U10,V10,XTIME,XLONG,XLAT ../D03_FLORENCE/WRF/wrfout_d01_2018-12-29_00:00:00.nc TEMPORARY/CTL_wrfout_d01_2018-12-29_00:00:00.nc
  NC_GLOBAL#HYBRID_OPT=2
  NC_GLOBAL#HYPSOMETRIC_OPT=2
  NC_GLOBAL#ICLOUD=1
  NC_GLOBAL#ICLOUD_CU=0
  NC_GLOBAL#IDEAL_CASE=0
  NC_GLOBAL#ISFFLX=1
  NC_GLOBAL#ISFTCFLX=0
  NC_GLOBAL#ISHALLOW=0
  NC_GLOBAL#ISICE=15
  NC_GLOBAL#ISLAKE=21
  NC_GLOBAL#ISOILWATER=14
  NC_GLOBAL#ISURBAN=13
  NC_GLOBAL#ISWATER=17
  NC_GLOBAL#I_PARENT_START=1
  NC_GLOBAL#JULDAY=363
  NC_GLOBAL#JULYR=2018
  NC_GLOBAL#J_PARENT_START=1
  NC_GLOBAL#KHDIF=0
  NC_GLOBAL#KM_OPT=4
  NC_GLOBAL#KVDIF=0
  NC_GLOBAL#MAP_PROJ=1
  NC_GLOBAL#MAP_PROJ_CHAR=Lambert Conformal
  NC_GLOBAL#MFSHCONV=0
  NC_GLOBAL#MMINLU=MODIFIED_IGBP_MODIS_NOAH
  NC_GLOBAL#MOAD_CEN_LAT=44.137997
  NC_GLOBAL#MOIST_ADV_OPT=2
  NC_GLOBAL#MP_PHYSICS=10
  NC_GLOBAL#NCO=4.4.4
  NC_GLOBAL#NTASKS_TOTAL=96
  NC_GLOBAL#NTASKS_X=8
  NC_GLOBAL#NTASKS_Y=12
  NC_GLOBAL#NUM_LAND_CAT=33
  NC_GLOBAL#OBS_NUDGE_OPT=0
  NC_GLOBAL#OPT_ALB=2
  NC_GLOBAL#OPT_BTR=1
  NC_GLOBAL#OPT_CROP=0
  NC_GLOBAL#OPT_CRS=1
  NC_GLOBAL#OPT_FRZ=1
  NC_GLOBAL#OPT_GLA=1
  NC_GLOBAL#OPT_INF=1
  NC_GLOBAL#OPT_PEDO=1
  NC_GLOBAL#OPT_RAD=3
  NC_GLOBAL#OPT_RSF=1
  NC_GLOBAL#OPT_RUN=3
  NC_GLOBAL#OPT_SFC=1
  NC_GLOBAL#OPT_SNF=4
  NC_GLOBAL#OPT_SOIL=1
  NC_GLOBAL#OPT_STC=3
  NC_GLOBAL#OPT_TBOT=1
  NC_GLOBAL#PARENT_GRID_RATIO=1
  NC_GLOBAL#PARENT_ID=0
  NC_GLOBAL#POLE_LAT=90
  NC_GLOBAL#POLE_LON=0
  NC_GLOBAL#PREC_ACC_DT=0
  NC_GLOBAL#RADT=2
  NC_GLOBAL#RA_LW_PHYSICS=4
  NC_GLOBAL#RA_SW_PHYSICS=4
  NC_GLOBAL#SCALAR_ADV_OPT=2
  NC_GLOBAL#SCALAR_PBLMIX=0
  NC_GLOBAL#SF_LAKE_PHYSICS=0
  NC_GLOBAL#SF_OCEAN_PHYSICS=0
  NC_GLOBAL#SF_SFCLAY_PHYSICS=2
  NC_GLOBAL#SF_SURFACE_MOSAIC=0
  NC_GLOBAL#SF_SURFACE_PHYSICS=4
  NC_GLOBAL#SF_URBAN_PHYSICS=2
  NC_GLOBAL#SGFDDA_END_H=0
  NC_GLOBAL#SGFDDA_INTERVAL_M=0
  NC_GLOBAL#SHCU_PHYSICS=0
  NC_GLOBAL#SIMULATION_INITIALIZATION_TYPE=REAL-DATA CASE
  NC_GLOBAL#SIMULATION_START_DATE=2018-12-29_00:00:00
  NC_GLOBAL#SKEBS_ON=0
  NC_GLOBAL#SMOOTH_OPTION=0
  NC_GLOBAL#SOUTH-NORTH_GRID_DIMENSION=141
  NC_GLOBAL#SOUTH-NORTH_PATCH_END_STAG=141
  NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG=140
  NC_GLOBAL#SOUTH-NORTH_PATCH_START_STAG=1
  NC_GLOBAL#SOUTH-NORTH_PATCH_START_UNSTAG=1
  NC_GLOBAL#SPEC_BDY_FINAL_MU=1
  NC_GLOBAL#SST_UPDATE=1
  NC_GLOBAL#STAND_LON=11.092
  NC_GLOBAL#START_DATE=2018-12-29_00:00:00
  NC_GLOBAL#SURFACE_INPUT_SOURCE=1
  NC_GLOBAL#SWINT_OPT=0
  NC_GLOBAL#SWRAD_SCAT=1
  NC_GLOBAL#TITLE= OUTPUT FROM *             PROGRAM:WRF-Chem V4.2.2 MODEL
  NC_GLOBAL#TKE_ADV_OPT=2
  NC_GLOBAL#TRACER_PBLMIX=1
  NC_GLOBAL#TRUELAT1=44.138
  NC_GLOBAL#TRUELAT2=44.138
  NC_GLOBAL#USE_Q_DIABATIC=0
  NC_GLOBAL#USE_THETA_M=0
  NC_GLOBAL#WEST-EAST_GRID_DIMENSION=146
  NC_GLOBAL#WEST-EAST_PATCH_END_STAG=146
  NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG=145
  NC_GLOBAL#WEST-EAST_PATCH_START_STAG=1
  NC_GLOBAL#WEST-EAST_PATCH_START_UNSTAG=1
  NC_GLOBAL#W_DAMPING=1
  NC_GLOBAL#YSU_TOPDOWN_PBLMIX=0
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":XLONG
  SUBDATASET_1_DESC=[140x145] longitude (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":XLAT
  SUBDATASET_2_DESC=[140x145] latitude (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":MEBIO_ISOP
  SUBDATASET_3_DESC=[8833x140x145] MEBIO_ISOP (32-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":PM10
  SUBDATASET_4_DESC=[8833x1x140x145] PM10 (32-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":PM2_5_DRY
  SUBDATASET_5_DESC=[8833x1x140x145] PM2_5_DRY (32-bit floating-point)
  SUBDATASET_6_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":SWDOWN
  SUBDATASET_6_DESC=[8833x140x145] SWDOWN (32-bit floating-point)
  SUBDATASET_7_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":T2
  SUBDATASET_7_DESC=[8833x140x145] T2 (32-bit floating-point)
  SUBDATASET_8_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":U10
  SUBDATASET_8_DESC=[8833x140x145] U10 (32-bit floating-point)
  SUBDATASET_9_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":V10
  SUBDATASET_9_DESC=[8833x140x145] V10 (32-bit floating-point)
  SUBDATASET_10_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":dvel_o3
  SUBDATASET_10_DESC=[8833x1x140x145] dvel_o3 (32-bit floating-point)
  SUBDATASET_11_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":no2
  SUBDATASET_11_DESC=[8833x1x140x145] no2 (32-bit floating-point)
  SUBDATASET_12_NAME=NETCDF:"wrfout_d03_Florence_Control.nc":o3
  SUBDATASET_12_DESC=[8833x1x140x145] o3 (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

As you can see, the variables have a 2d array for each spatial coordinate, and a 1d array for the time. Here is some code to see this:

import xarray as XR
ncfile = xr.open_dataset(file_path)
ncfile['T2']
Out[4]: 
<xarray.DataArray 'T2' (XTIME: 8833, south_north: 140, west_east: 145)>
[179309900 values with dtype=float32]
Coordinates:
  * XTIME    (XTIME) datetime64[ns] 2018-12-29 ... 2020-01-01
    XLONG    (south_north, west_east) float32 ...
    XLAT     (south_north, west_east) float32 ...
Dimensions without coordinates: south_north, west_east
Attributes:
    units:        K
    FieldType:    104
    MemoryOrder:  XY 
    description:  TEMP at 2 M

ncfile['T2']['west_east']
Out[7]: 
<xarray.DataArray 'west_east' (west_east: 145)>
array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
        14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
        28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
        42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
        56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
        70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
        84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
        98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
       112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
       126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139,
       140, 141, 142, 143, 144])
Dimensions without coordinates: west_east

ncfile['T2']['south_north']
Out[8]: 
<xarray.DataArray 'south_north' (south_north: 140)>
array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
        14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
        28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
        42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
        56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
        70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
        84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
        98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
       112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
       126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139])
Dimensions without coordinates: south_north

west_east and south_north are basically the indexes to reach the real coordinates, stored in the XLAT and XLONG variables: the value T2[t,i,j] is located at position XLAT[i,j], XLONG[i,j] at time XTIME[t]. So the number of spatial points, for each time step, is 145x140=20300.

What I would like to do is to take one time step figure of variable T2 and convert to GeoTIFF, using a code like:

    db = xr.DataArray(data=t0.values,
                  coords={"lat": (["x","y"], ncfile['XLAT'].values),
                          "lon": (["x","y"], ncfile['XLONG'].values)},
                  dims=["x","y"])
    db = db.rio.set_spatial_dims('y','x')

    db.rio.set_crs("epsg:4326")
    db.rio.to_raster(r"GeoTIFF.tif")

but I get this warning

[...]miniconda3/envs/spyder-geo/lib/python3.9/site-packages/rasterio/__init__.py:230: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific.
  s = writer(path, mode, driver=driver,

and the output Tiff file has no coordinates, as you can see here

Driver: GTiff/GeoTIFF
Files: GeoTIFF.tif
Size is 140, 145
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Lower Left  (       0.000,     145.000) (  0d 0' 0.01"E,145d 0' 0.00"N)
Upper Right ( 140.0000000,   0.0000000) (140d 0' 0.00"E,  0d 0' 0.01"N)
Lower Right (     140.000,     145.000) (140d 0' 0.00"E,145d 0' 0.00"N)
Center      (  70.0000000,  72.5000000) ( 70d 0' 0.00"E, 72d30' 0.00"N)
Band 1 Block=140x14 Type=Float32, ColorInterp=Gray

the dataArray db has the following properties

<xarray.DataArray (y: 145, x: 140)>
array([[284.9859 , 284.90915, 284.83188, ..., 275.89563, 275.94363,
        275.98633],
       [284.8853 , 284.8054 , 284.7248 , ..., 275.9049 , 275.95135,
        275.99228],
       [284.7894 , 284.70612, 284.62192, ..., 275.91513, 275.96002,
        275.99908],
       ...,
       [277.46445, 277.41336, 277.36038, ..., 275.68375, 275.73285,
        275.7803 ],
       [277.4038 , 277.35394, 277.30228, ..., 275.8227 , 275.8771 ,
        275.93002],
       [277.34274, 277.29395, 277.24344, ..., 275.95752, 276.01672,
        276.07465]], dtype=float32)
Coordinates:
    lat      (y, x) float32 43.26 43.27 43.28 43.29 ... 44.48 44.49 44.5 44.51
    lon      (y, x) float32 10.39 10.39 10.39 10.39 ... 12.19 12.19 12.19 12.19
Dimensions without coordinates: y, x

Answer 1

Finally solved with the following script

import netCDF4 as nc
import numpy as np
import rasterio
from rasterio.transform import from_origin

# domain of interest bounds
xmin = 11.17837
xmax = 11.3420135404019
ymax = 43.81461
ymin = 43.7324703002022

# Input NetCDF file path
file_path = "wrfout.nc"

# Reading file and extraxting coordinates
data = nc.Dataset(file_path)
latitudes = data.variables["XLAT"][:,:]
longitudes = data.variables["XLONG"][:,:]
# take a single time step
time0_T2 = data.variables["T2"][0,:,:]

# defining lists for the output dataset
lat, lon, temp, nc = [], [], [], []
# counter variable for cell indexing
count = 1
for i in range(latitudes.shape[0]):
    for j in range(latitudes.shape[1]):
        if ymin <= latitudes[i, j] <= ymax and xmin <= longitudes[i, j] <= xmax:
            lat.append(latitudes[i, j])
            lon.append(longitudes[i, j])
            # convert temperature from Kelvin to Celsius
            temp.append(time0_T2[i, j] - 273.15)
            nc.append(count)
            count = count + 1

# converting to numpy column array
A = np.transpose(np.array([nc,temp,lon,lat]))

# writing data to csv
np.savetxt('time0_T2.csv',A,delimiter=",",header="cell_number,T2_degC,LON,LAT")

# x and y lenghts. Need to round because the grid is not exactly vertical and horizontal
# but it is a little bit shifted. In the area of interest, it was manually checked that
# nx = 14 and ny = 10.
nx = np.unique(np.round(A[:,2],2)).shape[0]
ny = np.unique(np.round(A[:,3],3)).shape[0]
# x and y resolutions
deltax = (np.max(A[:,2]) - np.min(A[:,2]))/(nx-1)
deltay = (np.max(A[:,3]) - np.min(A[:,3]))/(ny-1)

# up left corner
minx = np.min(A[:,2])
maxy = np.max(A[:,3])

# set the transformation for rasterio. Note that the shift operated is in order to
# make the lat,lon points to be the center of each grid cell
transform = from_origin(minx-0.5*deltax,maxy+0.5*deltay, deltax, deltay)

# write it to tiff format
with rasterio.open('time0_T2.tiff', 'w', driver='GTiff', width=nx, height=ny, count=1, dtype=A[:,1].dtype, crs='EPSG:4326', transform=transform) as dst:
    dst.write(A[:,1].reshape(1, ny,nx))