Create a GeoTiff from CMAQ

CMAQ
1

I was recently asked how to make a GeoTiff from CMAQ data. I want the response to be searchable by others in the future, so I am posting it here.

You can convert NetCDF files to GeoTiff using gdal_translate.[1] The example below uses a 12US1 GRIDCRO2D to make a GeoTiff. As a test, I read it with QGis and add US Census boundary file and check that the overlay aligns correctly (it does see figure).

$ gdal_translate -a_ullr 0 299 459 0 -a_srs "+proj=lcc +lat_1=33.0 +lat_2=45.0 +lat_0=40.0 +lon_0=-97.0 +y_0=1728000.0 +x_0=2556000.0 +a=6370000 +b=6370000 +to_meter=12000.0 +no_defs" NETCDF:TestData/MCIP/GRIDCRO2D.12US1.35L.160101:PURB GRIDCRO2D_12US1_35L_PURB.tif

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute

Input file size is 459, 299

0...10...20...30...40...50...60...70...80...90...100 - done.

Notes

  • The a_srs parameter will depend on the domain you are using and can be figured out from I/O API metadata
    • lat_1=P_ALP; lat_2=P_BET; lon_0=XCENT; lat_0=YCENT;
    • y_0=-YORIG; x_0=-XORIG; to_meter=XCELL
    • a and b are the spheroid of the earth (usually the same);
  • The a_ullr and must be set to 0 NROWS NCOLS 0. Otherwise, the image will be upside down.
  • The NetCDF variable is referenced using NETCDF::

Below is an image using the created GeoTiff in QGis with a US Census State Boundary boundary to demonstrate that the result make sense.

image
image723×477 74.5 KB

I hope this is useful.

[1] gdal_translate — GDAL documentation

8 Likes
2

If you’re using a file with multiple time steps and layers, you’ll need to use -b to specify which band your using (e.g., -b 1). The appropriate number is an unraveling of the TSTEP and LAY dimensions. For example, if you’re file has 25 times and 35 layers the table below shows the progression

band time layer
1 1 1
2 1 2
3 1 3
… … …
36 2 1
… … …
400 12 15
… … …
873 25 33
874 25 34
875 25 35

You can use gdalinfo to get a full listing.

1 Like
3

Thank you! This is very helpful.