NetCDF and Python

Corvus

A version of this PyQGIS code was used to loop through NOAA time series data for The Great Lakes comprised of 28 NetCDF files with 365 layers in each one. The data was cropped and projected, then each individual layer had a custom colour ramp applied to it using a QGIS QML file before being output as a PNG image file. These became individual frames in a video with additional information added using both Inkscape and Blender.

Thanks to Underdark at AnitaGraser.com for providing the head-start for a Python newbie with PyQGIS101

(Note: this was first posted in 2023 and has been reposted, in part, to provide test content for this new blogging CMS)

#Loop through years
ice_yr = 'PATH_TO_NetCDF_file\\FILE_NAME.nc'
i_count = range(1995,2022,1)
for i_cnt in i_count:
    loop = "{}".format(i_cnt)
    ice_input = ice_yr[:95] + loop + '_FILE_NAME.nc'
    gerent = processing.run("gdal:warpreproject", {'INPUT':ice_input,'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:3857'),'RESAMPLING':0,'NODATA':None,'TARGET_RESOLUTION':None,'OPTIONS':'','DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':False,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']
    cliptif = processing.run("gdal:cliprasterbymasklayer", {'INPUT':gerent,'MASK':'PATH_TO_CLIP_POLYGON\\GEOPACKAGE.gpkg|layername=CLIPLAYER','SOURCE_CRS':None,'TARGET_CRS':None,'NODATA':None,'ALPHA_BAND':False,'CROP_TO_CUTLINE':True,'KEEP_RESOLUTION':True,'SET_RESOLUTION':False,'X_RESOLUTION':None,'Y_RESOLUTION':None,'MULTITHREADING':False,'OPTIONS':'','DATA_TYPE':0,'EXTRA':'','OUTPUT':'TEMPORARY_OUTPUT'})
    Fileout = 'PATH_TO_OUTPUT\\1995_1.png'
    count = range(1,366,1)
    for cnt in count:
        ist = "{}".format(cnt)
        bandx = processing.run("gdal:rearrange_bands", {'INPUT':cliptif['OUTPUT'],'BANDS':[cnt],'OPTIONS':'','DATA_TYPE':0,'OUTPUT':'TEMPORARY_OUTPUT'})
        styled = processing.run("native:setlayerstyle", {'INPUT':bandx['OUTPUT'], 'STYLE':'PATH_TO_QGIS_QML/ColourRamp.qml', 'OUTPUT':'TEMPORARY_OUTPUT'})
        QgsProject.instance().addMapLayer(styled['OUTPUT'])
        manager = QgsProject.instance().layoutManager()
        layout = manager.layoutByName("Export")
        exporter = QgsLayoutExporter(layout)
        Fileout = Fileout[:102] + loop + '_' + ist + '.png'
        print(Fileout)
        exporter.exportToImage(Fileout, QgsLayoutExporter.ImageExportSettings())
        project = QgsProject.instance()
        delete_this_layer = project.mapLayersByName('OUTPUT')[0]
        project.removeMapLayer(delete_this_layer.id())
        cnt = cnt + 1
    i_cnt = i_cnt + 1