Converting NetCDF to GeoTiff file using Python

Introduction

Here, we are going to use the Arctic sea ice concentration (SIC) data from the ASI (the Arctic Radiation And Turbulence Interaction Study (ARTIST) Sea Ice) at 6.25 km grid resolution in the Arctic. The ASI algorithm computes daily sea ice concentration from the polarization difference of 89 GHz brightness temperature. We are going to convert netCDF file to GeoTIFF using Python. GeoTIFF is based on the TIFF format and is used as an interchange format for georeferenced raster imagery.

Data

The ASI SIC data is available in netCDF format (.nc). Daily ASI SIC product is published by the University of Bremen. The ASI SIC product has a polar stereographic projection at the 6.25 km spatial resolution and referenced to the WGS-84 ellipsoid, with each grid point center latitude and longitude provided with the data. We will be using daily data of 1 February, 2020. The ASI SIC NetCDF file can be found here as “asi_AMSR2_n6250_20200201_v5.4.nc”.

Necessary Libraries

NetCDF files can be read with a few different Python modules. We will use netCDF4. Follow this article for more about using netCDF library. We will also use gdal, numpy, and maplotliblibraries. A little tip for installing gdal — create a new environment and install geopandas library using conda install geopandas and this will automatically install gdal as it is one of the dependencies. The above installation works for python 3.8.* versions.

Reading netCDF file

We will use the netCDF.Dataset() attribute for reading the netCDF file. Let’s read the netCDF file and explore the variables:

#Importing the necessary libraries
from osgeo import gdal
import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

#Reading the netCDF file
data = Dataset(r'\path\to\file.nc')
print (data)

The resultant output is:

Metadata of the netCDF file

Here, ‘x( =1216)’, ‘y (=1792)’ and ‘z (= 0 to 100)’ variables are the number of columns, rows and the values of SIC, respectively.

First, let’s try to display the netCDF file and see what it looks like just out of curiosity.

#Creating a vairable 'sic' to store the 2-D array variable 'z'
sic = data.variables['z'][:]

#Flipping the image upright in the axis = 0 i.e., vertically
asi_sic = np.flip(sic,0)

#Displaying the variable 'z' i.e., SIC
cmmap = plt.cm.get_cmap("jet")
cmmap.set_bad('dimgrey',1.)
fig, ax = plt.subplots()
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.05)

im = ax.imshow(asi_sic, cmap=cmmap, vmin=0, vmax=100)

fig.colorbar(im, cax=cax, orientation='vertical', label='Sea Ice Concentration (%)')

ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.gcf().set_size_inches(15, 15)
ax.set_title('ASI Daily SIC 1 February 2020 Grid 6.25 km', fontsize=20)
fig.set_dpi(300.0)
plt.show()

The resultant image is:

ASI SIC 1 February 2020

The reason why we flipped the sic variable is because the sic array is an inverted array. You may refer to this article for more details about this particular file.

Getting information for Georeferencing

First, we need a reference GeoTIFF file for georeferencing the netCDF file. So, we are going to use a GeoTiff file of the same netCDF file which is exported using QGIS software. We will read the georeferencing information using gdal as follows:

#Reading GeoTIFF file for GeoReferencing
ori_asi6 = gdal.Open (r"\path\to\file.tif")
gt_asi = ori_asi6.GetGeoTransform ()
proj_asi = ori_asi6.GetProjection()

We can see the variables as:

Geo Transform of the GeoTIFF file

Projection information of the GeoTIFF file

The gt_asi contains georeferencing information in the form of the six rows.

The second row is the X-direction pixel size (here in meters). Third and fifth rows are the rotational components which in this case are set to zero as the image is not rotated. The sixth row is the Y-direction pixel resolution (which is almost always negative). The first and fourth rows are the coordinates (here in meters) of the upper left corner pixel for the image.

The proj_asi variable contains the projection information which in our case is NSIDC Sea Ice Polar Stereographic North and the EPSG code is “3413” which is written in the last line at the end of the sentence.

Creating GeoTIFF file

We are now only left with creating the GeoTIFF file and is done using gdalas follows:

#Creating output file
driver = gdal.GetDriverByName("GTiff")#Getting GeoTIFF driver
driver.Register()

outds_asi = driver.Create ("ASI_625km_SIC_1_Feb_2020.tif", #Name of the putput file
                           xsize = asi_sic.shape[1], #Setting the number of columns
                           ysize = asi_sic.shape[0], #Setting the number of rows
                           bands = 1, #Setting the number of bands
                           eType = gdal.GDT_Float32) #Setting the data type i.e., float 32

outds_asi.SetGeoTransform(gt_asi) #Setting the geoinformation
outds_asi.SetProjection(proj_asi) #Setting the projection information
outband_asi = outds_asi.GetRasterBand(1) #Setting the number of bands
outband_asi.WriteArray(asi_sic) #Writing the 2-D asi_sic array on the file
outband_asi = None #Closing the file
outds_asi = None #Closing the file

Here, the output file is created using gdal driver. The output GeoTIFF file is named “ASI_625km_SIC_1_Feb_2020.tif” with the column and row numbers derived from the 2-D asi_sic array and data type as float 32. The geoinformation and the projection information are also set using gt_asi and proj_asi, respectively. The asi_sic 2-D array is written on the GeoTIFF file using WriteArray option. Finally, closing the files by setting outband_asi and outds_asi to None. This is crucial as the data will not be properly written on the output file without closing the file. And, that’s it!!!

Let’s check if the output GeoTIFF file is properly created or not by opening the file using QGIS.

Raster information of the GeoTIFF file created.

As we can see in the raster information and CRS sections that the information that we provided are correctly written on the file.

Conclusion

We have just learnt how to convert netCDF file to GeoTIFF file using python. GeoTIFF file is much easier to deal with as it contains spatial referencing information in the form of georeferenced raster imagery. Of course, we can use GIS softwares like QGIS, ArcGIS, etc. for creating GeoTIFF file from netCDF but using python, we get to know in depth how the whole process is being done. This is the key tenet behind this article.