用R语言计算GRIB格式气候数据的月均值

hillside

Compute Monthly Averages for GRIB Format Climate Data with R
http://www.nceas.ucsb.edu/scicom ... teGRIBAveragesWithR

This case demonstrates:
Extract monthly climate parameter estimates from General Circulation Model (GCM) model results from GRIB-format data files.
Aggregate the monthly data into long-term monthly average values using an R script
Write the output data into ArcMap GIS-compatible ASCII GRID format files
Import and display the data usingR, MATLAB and ArcMap GIS.
Download data and R Code for this example

Scenario: An NCEAS associate required global-scale estimates of three climate parameters (precipitation, humidity, and temperature) for two 10-year time periods 1981 - 1990 and 2081 - 2090. One good source for such data are Coupled General Circulation Models (CGCM) such as the HADCM3 model developed and managed by the United Kingdom Meteorological Office Hadley Centre for Climate Change.

After some investigation, the associate identified a suitable data set, generated using the Hadley Centre's HADCM3 model and distributed by The IPCC Data Distribution Centre using The GRIdded Binary or GRIB file format created by the World Meteorological Organization (WMO). GRIB is is a machine-independant binary data format optimized to support fast computer-to-computer transmission of large volumes of raster (grid)-format data such as digital Earth images or large scale geophysical simulation model results. In the current example. Output from HADCM3 General Circulation Model typically includes multiple grid layers for long time series (as many as 100 years) at a monthly temporal resolution, covering the entire Earth using a 72 row, 96 column grid (definition) with a resolution of 3.75 degrees (longitude) by 2.5 degrees (latitude).

The associate decided to 1) acquire these data sets for the desired parameters, 2) aggregate each parameter's monthly estimate into long-term monthly averages (e.g., average June precipitation for 1981 - 2081) and 3) write the final results into the ASCII ArcMap GRID raster grid format. The processing was performed using the open-source wgrib utility available from the National Weather Service, and an R Programming Language script. The associate wished to initially display the data in ArcMap GIS, and then export them into MATLAB for further analysis.

We demonstrate the techniques using the mean monthly temperature parameter.

Step 1: Discover and Download GRIB format GCM data sets.
We will use IPCC Data Distribution Center download site to identify and download relative humidity time series data. For this example, we select data from the SRES A2 scenario produced with HADCM3,and the time periods 1981 - 1990 and 2081 - 2090.

The download site delivers the GRIB-format data in compressed 'gzip' / Unix 'tar' archive format. Unix users use the unzip utility to extract the archive; MS Windows platform can use the 7-zip archive utility to do so. As indicated above, the output for each climate parameter from HADCM3 is a two dimensional grid covering the entire Earth surface and containing 96 columns and 72 rows. Each GRIB file contains the data for one climate parameter and all of the months for one 10-year time period: Each file therefore contains 120 raster 'layers', one for each month.

Step 2: Extract desired climate model parameters using the wgribArc.c utility
Next, we extract the data from the GRIB file into two-dimensional ASCII raster data files.

Several GRIB file data extraction software packages are available to researchers, from research centers such as the University Corporation for Atmospheric Research (UCAR), and The Institute of Global Environment and Society, Inc. (IGES), and the British Atmospheric Data Centre (BADC).For MATLAB users the OpenML routine read_grib that integrates with MATLAB, allows MATLAB scripts to decode GRIB files. read_grib is a MATLAB script 'wrapper' for the wgrib utility.

This project requires that the data in the GRIB file be written in ArcMap GIS-compatible ASCII GRID format. Unfortunately, there are no open source GRIB file utilities that write ArcMap GRID-format files (and in fact only one commercial package). Thus, we needed to modify the wgrib C language source code to add the six-line header:

NCOLS 96

NROWS 72

XLLCORNER -180.000000

YLLCORNER -90.000000

CELLSIZE 3.125000

NODATA_VALUE 0.000000

Call this modified version 'wgribArc'. Click here to review the wgribArc.c source code. Search for the string ' // NCEAS change' to review the modifications made.

We execute wgribArc one time time to for each GRIB file we wish to process. Command line parameters control the extraction process. Typically, we use a batch file to contain the wgrib command line, which can be lengthy. Here is a summary of the command line parameters used in each line of the batch file:

-g <filename>: path name to input GRIB file

-d ALL: Extract ALL records from time climate simulation time series

-yf 1981: First year of time series to extract

-yl 1990: Last year of time series to extract

-text: Output file will be a TEXT file

-S A2: Suffix to be appended to (automatically-generated) output file 'root name' .

-R: Output file root path name

-vt: type of output parameter that is being extracted (e.g., temp,wind)

Click here to view the sample batch file used to extract temperature data from two GRIB files.

The batch file calls wgribArc twice, once for each GRIB file (each covers one 10-yr time period), and produces 240 separate ASCII grid files, one for each month in each 10-year time period. These data files will be consolidated into 24 long-term monthly average files in Step 3.

Step 3: Create long-term monthly average temperature grids using R
Next, we generate a single 10-year average temptress grid for each month and each of the two 10-year time periods.

Write an R script that:

1) Computes a long-term average monthly temperature grid using the 10 monthly data files for each month

2) Writes the long-term average grid to an ArcMap ASCII GRID file,

3) Produces (and writes to a PDF file) a set of histograms for each monthly temperature file as well as the 10-year average file.

Click here to display the file masterfile.txt

Click here to display one of the input files generated by wgribArc, listed in masterfile.txt

Click here to display one of the individual monthly average files (June, 1990).

Click here to display one of the long-term monthly average files (June,1981 - 1990).

Click here to view the diagnostic plots for the 10 June averages,years 1981 - 1990.

We present two versions of this script: ProcGridAscFile.r using standard R data frames, and ProcGridAscFileMapTools.r, which uses the R maptools package methods for reading ArcMap GRID files.

To execute this step, start the R environment, then from the command prompt:

> library(maptools)

> source("ProcGridAscFileMapTools.r")

> MakeLongTermAvgClimateGrids(nFiles) # nFiles: number of lines in GRIB batch file (Step 1).

Step 4: Display the ASCII raster files using R, MATLAB, and ArcMap GIS
Using base R graphics functions:
View output

#
# Read the ASCII Grid into a SpatialGridDataFrame object
# Assume that current working directory has the following file:
#
WindImage = readAsciiGrid("wind_198190av06.asc")
#
# create a 32-level grey scale
#
breaks = seq(0,1,0.03125) # 32 levels between 0 and 1
        scale = grey(breaks)
#
# display the image with the scale
#
image(WindImage,col = scale)
#
Using MATLAB Mapping Toolbox:
View output

        [Z,R] = arggridread('wind_1980av06.asc')
        mapshow(Z,R,'displayType','surface')
Using ArcMap GIS Toolbox: Import ASCII GRID

Display Climate Data Image

Learning More:
GRIB File Format
National Center for Atmospheric Research: A Guide to WMO GRIB

Comprehensive History of the GRIB File Format

GRIB File Extraction and Decoding Software
University Corporation for Atmospheric Research (UCAR): GRIB and NetCDF Decoders

University of North Carolina: OpenML read_grib decoder for MATLAB

Institute for Global Environment and Society (IGES): GrADS

British Atmospheric Data Centre (BADC): GRIBEX

General Circulation Model (GCM) Information and Data Sets
Intergovernmental Panel on Climate Change (IPCC) IPCC

IPCC Data Distribution Center Home Page

Lawrence Livermore National Laboratory Program for Climate Model Diagnosis and Intercomparison

World Climate Research Program/ Coupled Model Intercomparison Project (WCRP / CMIP3)

National Center for Atmospheric Research

Earth and Sun Systems Laboratory (ESSL) Climate and Forecasting Models

Point of Contact for this Use Case: Rick Reeves, NCEAS Scientific Programmer
This Use Case compiled May, 2007 / updated October,2009

scicomp:
Spatial