0

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

1 Answer 1

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))
1
  • Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.
    – Community Bot
    Commented Jun 7, 2023 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.