You have used Python with GDAL and NumPy in order to create a command-line utility to import a NetCDF dataset into PostGIS.

A NetCDF dataset is composed of multiple subdatasets, and each subdataset is composed of multiple raster bands. Each band is composed of cells. This structure should be clear to you after investigating a sample NetCDF dataset using the gdalinfo GDAL command tool.

There are several approaches to exporting cell values to PostGIS. The approach you adopted here is to generate a PostGIS point layer for each subdataset, which is composed of one field for each subdataset band. You then iterated the raster cells and appended a point to the PostGIS layer with the values read from each cell band.

The way you do this with Python is by using the GDAL Python bindings. For reading, you open the NetCDF dataset, and for updating, you open the PostGIS database, using the correct GDAL and OGR drivers. Then, you iterate the NetCDF subdatasets, using the GetSubDatasets method, and create a PostGIS table named NetCDF subdataset variable (with the prefix) for each subdataset, using the CreateLayer method.

For each subdataset, you iterate its bands, using the GetRasterBand method. To read each band, you run the ReadAsArray method which uses NumPy to get the band as an array.

For each band, you create a field in the PostGIS layer with the correct field data type that will be able to store the band's values. To choose the correct data type, you investigate the band's data type, using the DataType property.

Finally, you iterate the raster cells, by reading the correct x and y coordinates using the subdataset transform parameters, available via the GetGeoTransform method. For each cell, you create a point with the CreateGeometryFromWkt method, then set the field values, and read from the band array using the SetField feature method.

Finally, you append the new point to the PostGIS layer using the CreateFeature method.