3

I need to convert netCDF files from AWS's GOES data into GeoTIFF. My end goal is to use these as raster layers in GeoServer. I have used GDAL to warp a GeoTIFF into EPSG 4326 so for now I just need to get it into a GeoTIFF. I tried using gdal_translate with the command:

gdal_translate NETCDF:"OR_ABI-L2-CMIPF-M3C02_G16_s20180640000387_e20180640011153_c20180640011226.nc":CMI test2.tif

This gives me an image, but opening it with an image viewer just shows a white background with a completely black circle. I've also tried using Python with the code:

# Required libraries
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import pdb

# Path to the GOES-R simulated image file
path = '/home/brian/projects/aws_goes_testing/OR_ABI-L2-CMIPF-M3C02_G16_s20180640000387_e20180640011153_c20180640011226.nc'

# Open the file using the NetCDF4 library
nc = Dataset(path)

# Extract the Brightness Temperature values from the NetCDF
data = nc.variables['CMI'][:]

# Show data
plt.imshow(data)
plt.show()

This gives me an image I want but I'm not sure how I can save the data variable as a GeoTIFF. Obviously I can use plt.savefig() but this isn't what I want since it will include the axis and everything from the figure I don't want.

2

I managed to get it to kind of work. When I put the geotiff in geoserver and load it onto the cesium globe it doesn't line up properly, so something isn't write with the projection, but the general idea is there:

#!/usr/bin/python3

# Required libraries
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from osgeo import gdal
from osgeo import osr
import numpy as np
import pdb
import subprocess

# Path to the GOES-R simulated image file
path = '/home/brian/projects/aws_goes_testing/OR_ABI-L2-CMIPF-M3C02_G16_s20180661800402_e20180661811169_c20180661811241.nc'

# Open the file using the NetCDF4 library
nc = Dataset(path)

# Extract the Brightness Temperature values from the NetCDF
data = nc.variables['CMI'][:]
#it automatically fills null values to 65535
data.data[data.data == 65535] = -1
normData = data.data * 255
normData[normData == -255] = 255
normData8 = normData.astype('uint8')

nx = data.data.shape[0]
ny = data.data.shape[1]
latLon = nc['geospatial_lat_lon_extent']
xmin, ymin, xmax, ymax = [latLon.geospatial_westbound_longitude, latLon.geospatial_southbound_latitude, latLon.geospatial_eastbound_longitude, latLon.geospatial_northbound_latitude]


xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)
normData = data.data/data.data.max() * 255
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF_unwarped2.tif', ny, nx, 1, gdal.GDT_Byte)

dst_ds.SetGeoTransform(geotransform)    # specify coords
srs = osr.SpatialReference()            # establish encoding
srs.ImportFromEPSG(4326)                # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(normData8)   # write r-band to the raster
dst_ds.FlushCache()                     # write to disk

dst_ds = None
# Show data
#plt.imshow(data, cmap='Greys')
#plt.show()
nc.close()

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.