🧐"Not-so-good" NetCDF

Sometimes you can encounter problems with data georeferencing, as metadata could be written not in a QGIS-friendly style.

Data source

ESA Permafrost Climate Change Initiative, Permafrost extent for the Northern Hemisphere, v3.0

QGIS

Adding data to the map as usual: the layer has incorrect georeference.

data view in QGIS

Python + XARRAY + GDAL

Let's try to study the dataset with Xarray

The full notebook download is below. To run it, set up conda environment, savethe notebook to your working directory, and run Jupyter Notebook from this directory.

100KB

Permafrost_NETCDF_exp.ipynb

Read and view data

import xarray as xr

ds = xr.open_dataset('ESACCI-PERMAFROST-L4-PFR-MODISLST_CRYOGRID-AREA4_PP-2019-fv03.0.nc')

Dataset has 2 variables, one of which is crs, but x and y coordinates are presented as 1-d arrays

prf = ds.PFR[0,:,:]
plt.imshow(prf, cmap='Spectral', vmin=0, vmax=100)
plt.colorbar()

data view: the rotation is correct

To export data as GeoTIFF some info georeferencing is needed. For GDAL it is a projection (crs reference) and geotransform parameters (upper left x and y coordinates and pixel size). Some of this info is stored in the polar_stereographic variable.

crs_wkt is in GDAL-compatible format, we can use it. But GeoTransform is not consistent with dataset upper left pixel coordinates:

print(ds.polar_stereographic.GeoTransform)
print(ds.x.values.min(), ds.x.values[1] - ds.x.values[0], 0, ds.y.values.max(), 0, ds.y.values[1] - ds.y.values[0])
> -8679599.425969256 926.6254331383326 0 4120958.560533902 0 -926.6254331383326 
> -6111475.222394747 926.6254331385717 0 4114895.094696621 0 -926.625433138106

Save to GeoTIFF

def save_file_singleband(input_array, out_file, projection, geotransform):
    
    """Save numpy array to singleband GTiff file"""
    
    driver = gdal.GetDriverByName('GTiff')
    
    dataset = driver.Create(out_file,
                        input_array.shape[1],
                        input_array.shape[0],
                        eType=gdal.GDT_Float32)
    
    #set projection
    dataset.SetProjection(projection)
    # set geotransform
    dataset.SetGeoTransform(geotransform)

    outband = dataset.GetRasterBand(1)
    outband.WriteArray(input_array)
    
    dataset = None

save_file_singleband(
    prf,
    'permafrost_proj_wkt.tif',
    ds.polar_stereographic.crs_wkt,
    [-6111475.22239475,926.62543314,0,4114895.09469662,0,-926.62543314]
)

Output file in QGIS map view