1

I have netCDF data I convert to geoTiff using the following code:

import numpy as np
from osgeo import gdal
from osgeo import osr
from netCDF4 import Dataset
import glob


def open_netCDF(filename, path):

    """
    Opens netCDF and returns dataset
    """

    file = Dataset(path + filename + ".nc")

    # get variables
    time = file.variables["time"]
    lon = file.variables["lon"]
    lat = file.variables["lat"]
    data = file.variables["data"]

    # three variables as numpy arrays
    lonarray = np.array(lon)
    latarray = np.array(lat)
    dataarray = np.array(data)

    return dataarray, latarray, lonarray


def write_GeoTiff(data, lats, lons, filename):

    # get bounds
    xmin,ymin,xmax,ymax = [lons.min(),lats.min(),lons.max(),lats.max()]

    # no. rows and columns
    nrows,ncols = np.shape(data)

    # resolution of x/y
    xres = (xmax-xmin)/float(ncols)
    yres = (ymax-ymin)/float(nrows)

    geotransform=(xmin, xres, 0,ymax ,0, -yres)

    output_raster = gdal.GetDriverByName('GTiff').Create(filename + ".tiff",ncols, nrows, 1 ,
                                                         gdal.GDT_CInt16)  # Open the file

    output_raster.SetGeoTransform(geotransform)            # Specify its coordinates
    srs = osr.SpatialReference()                           # Establish its coordinate encoding
    srs.ImportFromEPSG(4326)                               # specifies WGS84 lat lon
    output_raster.SetProjection(srs.ExportToWkt())         # Exports coordinate system
                                                           # to the file
    output_raster.GetRasterBand(1).WriteArray(data)        # Writes my array to the raster
    output_raster.FlushCache()
    del output_raster

if __name__ == '__main__':

    ncfiles = glob.glob("/example/path/*.nc")

    for ncfile in ncfiles:

        ncfile = ncfile.split('/')
        ncfile = ncfile[-1]
        ncfile = ncfile.split('.')
        ncfile = ncfile[0]

        dataarray, latarray, lonarray = open_netCDF(ncfile, "example/path/")

        write_GeoTiff(dataarray, latarray, lonarray, ncfile)

My code runs without complaint. But, when I open in a image viewer (ImageMagick for example) it is black.

I can open with ArcGIS no problem, it works fine.

What is ArcGIS reading/defining that ImageMagick can't from the source GeoTiff?

5
  • 3
    In geospatial images with high bit depth the data is very often using just a narrow part of the whole range that is supported by 15, 32, or 64 bits. If all data is close to the black end of the supported range and viewer is showing the whole range the image appears Black to eye. Geospatial viewers usually check the min and max values of data from the histogram and stretch the data range to show as min and max lightness on screen. ImageMagick probably does not do this contrast enhancement automatically.
    – user30184
    Feb 20 '20 at 15:06
  • Hi - thanks for your response. I'm not sure that's what's going on here as; When I dump the pixel values to a text file all the lines display: ( 0, 0, 0) #000000000000 gray(0,0,0). So it seems like there is no colour data?? [to dump pixels: convert image.tiff pixels.txt]
    – dwright
    Feb 20 '20 at 16:35
  • 1
    GeoTiffs rarely open with none GIS software, I'm not sure if it is the extra metadata (origin, projection etc) or the compression or tiling that breaks them but it does.
    – Ian Turton
    Feb 20 '20 at 17:00
  • 1
    You can check the min/max values also with gdalinfo -stats image.tif. If ArcGIS opens the image it is unlikely that all pixel values are zeroes. Have you tried to open the image with QGIS?
    – user30184
    Feb 20 '20 at 17:18
  • Ah, I see the expected stat from gdalinfo -stats, but, the tiff appears blank in QGIS?!
    – dwright
    Feb 21 '20 at 11:33
1

I have recently used part of your code (all of it tbh), to finally finish my NetCDF to GeoTIFF script but I had problems with QGIS not displaying the data, so one of the changes I made was the way I get the type of the data I'm trying to display.

So I changed this, where you are manually setting the type of the data (QGIS wasn't reading GDT_CFloat32 for some reason):

output_raster = gdal.GetDriverByName('GTiff').Create(filename + ".tiff",ncols, nrows, 1 ,
                                                     gdal.GDT_CInt16)

To this (where the type was just Float32):

data_type_code = gdal_array.NumericTypeCodeToGDALTypeCode(data.dtype)
output_raster = gdal.GetDriverByName('GTiff').Create(
    filename + ".tiff", ncols, nrows, 1, data_type_code)

I've been able to open them using IrfanView on Windows, but before I made the previous change, it popped an error saying the TIFF data type was not supported.

I think that's the reason you can't display them, because ImageMagick does not support the data types embedded in the TIFF.

2
  • Hi - I actually posted another question trying to resolve why my .tiff would not open in QGIS. I followed your method and this worked perfectly (in QGIS) - thank you!
    – dwright
    Feb 24 '20 at 14:13
  • You should either post the answer there, or redirect it to here so other people who have the same problem can actually see a posible answer instead of an unanswered one.
    – Parrawk
    Feb 24 '20 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.