netcdf to .img or geotif format

Miscellaneous questions you have about anything related to graphically displaying WRF output

netcdf to .img or geotif format

Postby Bilal » Sat May 19, 2012 10:33 am

Hi all,

I want to convert netcdf (WRF output) file into .img or .geotif format.

Is anyone know how to convert it?

Thank you.

Regards,
Bilal
Bilal
 
Posts: 14
Joined: Sun Feb 05, 2012 1:47 pm

Re: netcdf to .img or geotif format

Postby nh_modeler » Tue Sep 06, 2016 6:12 pm

Here's an example of an R script I use to convert netcdf data to geotiff format. The example creates a geotiff of the HGT_M field from geogrid output. It works well for displaying domains with Lambert projection in a GIS framework.

Code: Select all
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# FILE: WRFnc_to_tif.R
# AUTHOR: Sandra LeGrand ERDC-CRREL-NH
# DATE CREATED: Spring 2016
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#
# PURPOSE
# =======
# WRFnc_to_tif.R converts WRF geogrid netcdf data in Lambert Conic
# Conformal projection to geotiff format.
#
# Output file uses naming convention VARNAME.tif.
#
# RUN INSTRUCTIONS
# ================
# source("~/WRFnc_to_tif.R")
#
# Note: User must set configuration vars below.
#
# oDir     -> Output directory path
# gDir     -> Geogrid data dirctory path
# filename -> Geogrid file name
# varname  -> User selected two-dimensional geogrid variable
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#----------------------------------------------------------------------
# Set hard-coded parameters.
#----------------------------------------------------------------------

oDir     <- '.'
gDir     <- '../DATA/STATIC/'
filename <- 'geo_em.d02.nc'
varname  <- 'HGT_M'

#----------------------------------------------------------------------
# Load required libraries.
#----------------------------------------------------------------------

library('ncdf')                 # NetCDF
library('rgdal')                # R Geospatial Data Abstraction
library('sp')                   # Spatial Data
library('raster')               # Raster I/O

#----------------------------------------------------------------------
# Open geogrid file.
#----------------------------------------------------------------------

# get file path
filename <- paste0(gDir, '/', filename)

# open geogrid file
ncid <- open.ncdf(filename)

#----------------------------------------------------------------------
# Set projection info; assume lambert conic conformal.
#----------------------------------------------------------------------

# pull global attributes for projection info
lon_0    <- att.get.ncdf(ncid,0, 'STAND_LON')[[2]]
lat_0    <- att.get.ncdf(ncid,0, 'MOAD_CEN_LAT')[[2]]
truelat1 <- att.get.ncdf(ncid,0, 'TRUELAT1')[[2]]
truelat2 <- att.get.ncdf(ncid,0, 'TRUELAT2')[[2]]


# generate projection parameter string (rgdal uses PROJ.4 parameters)
str <- paste0('+proj=lcc',           # Projection name
              ' +lat_1=', truelat1,  # Lat of 1st standard parallel
              ' +lat_2=', truelat2,  # Lat of 2nd standard parallel
              ' +lat_0=', lat_0,     # Outermost domain grid center lat
              ' +lon_0=', lon_0,     # Standard meridian
              ' +datum=WGS84',       # 1984 World Geodetic System Datum
              ' +ellps=sphere',      # WRF assumes globe is a sphere
              ' +a=6370000',         # Semimajor ellipsoid axis     
              ' +b=6370000',         # Semiminor ellipsoid axis
              ' +units=m',           # Units in meters, not degrees
              ' +no_defs')           # Ignore any default settings

# generate corner point projection string
str2 <- paste0('+proj=longlat',      # Projection name
              ' +a=6370',            # Semimajor ellipsoid axis
              ' +b=6370',            # Semiminor ellipsoid axis
              ' +towgs84=0,0,0',     # Datum transform parameters
              ' +no_defs')           # Ignore any default settings

# set the coordinate reference system arguments
crs.lcc <- CRS(str)
crs.ll <- CRS(str2)

#----------------------------------------------------------------------
# Set map extent info.
#----------------------------------------------------------------------

# pull WRF corner point arrays from global attributes
corner_lons <- att.get.ncdf(ncid,0, 'corner_lons')[[2]]
corner_lats <- att.get.ncdf(ncid,0, 'corner_lats')[[2]]

# pull corner lat/lons from unstaggered grid for plot extent
ll <- data.frame()
ll[1,1] <- corner_lons[13] # unstaggered grid lower left corner lon
ll[1,2] <- corner_lats[13] # unstaggered grid lower left corner lat
ll[2,1] <- corner_lons[15] # unstaggered grid upper right corner lon
ll[2,2] <- corner_lats[15] # unstaggered grid upper right corner lat

# create spatial points from coordinates data frame
ll <- SpatialPoints(ll, crs.ll)                  # must do this to flip
                                                 # degrees into meters
# transform them to special WRF projection
ll <- spTransform(ll, crs.lcc)

#----------------------------------------------------------------------
# Read in parameter.
#----------------------------------------------------------------------

# create dataframe
var <- data.frame(get.var.ncdf(ncid, varname))
varsize  = dim(var)

# reverse array
var <- var[,rev(colnames(var))]

# transpose array and turn into matrix
var <- t(data.matrix(var, rownames.force = NA))

#----------------------------------------------------------------------
# Convert matrix into a raster and apply georeferencing info.
#----------------------------------------------------------------------

x <- raster(var)
projection(x) <- crs.lcc
extent(x) <- extent(ll)

#----------------------------------------------------------------------
# Write out data in geotiff format.
#----------------------------------------------------------------------

# create output file name and path
out_file = paste0(oDir, '/', varname, '.tif')

# write output faster
writeRaster(x, filename=out_file, format='GTiff', overwrite=TRUE)
nh_modeler
 
Posts: 2
Joined: Thu Aug 08, 2013 7:05 pm


Return to Miscellaneous

Who is online

Users browsing this forum: No registered users and 2 guests