Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Extracting specific netcdf info and converting to GeoTIFF in python

I am attempting to extract a specific set of data from a netCDF file and then convert said data to a GeoTIFF.

So far I have managed to extract the data I want using netCDF4, all the data in the file are stored as 1d arrays (lat, lon, data I want) and assigning them to a 2d array. The netcdf file that I am working with was subsetted to a specific region. From here however I am at a loss.

I have a slight understanding of how geotiff conversion works via what I have read at these links:

https://borealperspectives.wordpress.com/2014/01/16/data-type-mapping-when-using-pythongdal-to-write-numpy-arrays-to-geotiff/

http://adventuresindevelopment.blogspot.co.uk/2008/12/create-geotiff-with-python-and-gdal.html

And here is what I have currently:

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

#Reading in data from files and extracting said data
ncfile = netCDF4.Dataset("data.nc", 'r') 
dataw = ncfile.variables["dataw"][:]
lat = ncfile.variables["Latitude"][:]
long = ncfile.variables["Longitude"][:]


n = len(dataw)
x = np.zeros((n,3), float)

x[:,0] = long[:]
x[:,1] = lat[:]
x[:,2] = dataw[:]

nx = len(long)
ny = len(lat)
xmin, ymin, xmax, ymax = [long.min(), lat.min(), long.max(), lat.max()]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)

#Creates 1 raster band file
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 1, gdal.GDT_Float32)

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

The process of geotiff creation above I have largely sourced from here:

How do I write/create a GeoTIFF RGB image file in python?

I have attempted this based on the understanding that the array I want to write to the raster is a 2d array with 3 columns, 2 of co-ords and 1 of data. My result which I check in snap is a black page with a white line along the LHS.

So my question is as follows:

How can I extract the necessary geotransform data from my netcdf file, adjust the geotransform parameters appropriately and subsequently use the extracted lat long + dataw arrays to write to a geotiff file?

like image 260
cd123 Avatar asked Nov 07 '22 23:11

cd123


1 Answers

Try using gdal_translate command line to convert to geotiff file

gdal_translate NETCDF:"<filename.nc>":<varible name> <required_file>.tif
like image 158
Shreya Gupta Avatar answered Nov 14 '22 21:11

Shreya Gupta