G5NR Data Access Guide: Difference between revisions

From GEOS-5
Jump to navigation Jump to search
Pchakrab (talk | contribs)
Pchakrab (talk | contribs)
Line 182: Line 182:
  import netCDF4 as nc4
  import netCDF4 as nc4
   
   
  rootgrp = nc4.Dataset('http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/in\
  rootgrp = nc4.Dataset('http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv', 'r')
st01hr_3d_T_Cv', 'r')
   
   
  # read air temperature                                                                                 
  # read air temperature                                                                                 

Revision as of 10:45, 16 October 2014


Background

File spec

Model config

Getting data

ftp, http

Download tool

opendap

Client access

In the following, we read the field 'T' (air temperature) from collection http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv.

Programming

These are simple programs to read the air temperature and compute its min/max. These codes require an opendap enabled NetCDF-4 library. The utility nc-config (nf-config for Fortran) bundled with the NetCDF-4 installation can be used to determine the necessary compiler flags.

C

#include<stdio.h>
#include<stdlib.h>
#include<netcdf.h>  // for reading NR files

/* Handle errors by printing an error message and exiting with a
 * non-zero status. */
#define ERRCODE 2
#define ERR(e) {printf("Error: %s\n", nc_strerror(e)); exit(ERRCODE);}

int main(void){
  // file name
  char* T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv";
  
  // netCDF ID for the file and data variable
  int ncid, varid;

  // global 4D array: (time,lev,lat,lon), one time step
  const int IM = 720;
  const int JM = 361;
  const int LM = 72;
  const int asyz = 1*LM*JM*IM;
  float *T = NULL;

  // hypercube for reading one array for a given time
  size_t start[4] = {36, 0, 0, 0}; // time step 37
  size_t count[4] = {1, LM, JM, IM}; // 1 time step, 3D (lon,lat,lev) array

  // return code
  int rc;

  // min/max values
  float minval, maxval;

  // misc counters
  int ctr;

  // allocate memory for T
  T = malloc(asyz*sizeof(float));

  // read the data file
  printf("Reading T..."); fflush(stdout);
  if (rc = nc_open(T_file, NC_NOWRITE, &ncid)) ERR(rc);
  if (rc = nc_inq_varid(ncid, "t", &varid)) ERR(rc);
  if (rc = nc_get_vara_float(ncid, varid, start, count, T)) ERR(rc);
  printf("done.\n"); fflush(stdout);

  // min/max of T
  minval = 1.0e15;
  maxval = -1.0e15;
  for (ctr=0; ctr<asyz; ctr++){
    if (T[ctr]<minval){
      minval = T[ctr];
    } else if (T[ctr]>maxval){
      maxval = T[ctr];
    }
  }
  printf("min(T): %f\n", minval);
  printf("max(T): %f\n", maxval);

  // free memory
  free(T);
  
  return 0;
}

Fortran

program g5nr_reader

  use netcdf           ! for reading the NR files                                                                               

  implicit none

  !  File name                                                                                                                 
  !  ---------                                                                                                                 
  character(len=256) :: T_file

  !  Global, 4D array: (lon,lat,lev,time)                                                                                      
  !  ------------------------------------ 
  real, pointer :: T(:,:,:,:) => null()

  !  Miscellaneous                                                                                                              
  !  -------------                                                                                                              
  integer :: ierr
  integer :: im, jm, lm
  integer :: ncid, varid
  integer :: start(4), count(4)

  !  For now hard code file name and dimensions                                                                                
  !  ------------------------------------------                                                                                
  im = 720
  jm = 361
  lm = 72
  T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv"

  !  Allocate the Global 4-D array with only 1 time level                                                                      
  !  ----------------------------------------------------                                                                      
  allocate(T(im,jm,lm,1))

  !  Hypercube for reading one 1 array for a given time                                                                         
  !  -------------------------------------------------                                                                          
  start = (/  1,  1,  1, 37 /)   ! time level 37                                                                                
  count = (/ im, jm, lm, 1  /)   ! 1 time level, 3D (lon,lat,lev) array                                                         

  !  Read the data file                                                                                                         
  !  ------------------                                                                                                         
  write(*,*)'Reading T'
  call check( nf90_open(T_file,NF90_NOWRITE,ncid), "opening T file")
  call check( nf90_inq_varid(ncid,"t",varid), "getting T varid")
  call check( nf90_get_var(ncid,varid,T,start=start,count=count), "reading T")
  call check( nf90_close(ncid), "closing T file")

  !  Orint min/max of arrays                                                                                                    
  !  -----------------------                                                                                                    
  write(*,*)'T: ', maxval(T),minval(T)

  !  All done                                                                                                                   
  !  --------                                                                                                                   

contains

  subroutine check(status, loc)

    integer, intent(in) :: status
    character(len=*), intent(in) :: loc

    if(status /= NF90_NOERR) then
       write (*,*) "Error at ", loc
       write (*,*) NF90_STRERROR(status)
    end if

  end subroutine check

end program g5nr_reader

Shmem example

Free clients

Python
netcdf4-python

If netcdf4-python module is available, the following script would read and compute the min/max value of air temperature for the specified time.

#!/usr/bin/env python                                                                                 

import sys
import numpy as np
import netCDF4 as nc4

rootgrp = nc4.Dataset('http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv', 'r')

# read air temperature                                                                                
print 'Reading T for time=37...',; sys.stdout.flush()
Ttime37 = rootgrp.variables['t'][36,:,:,:]
print 'done.'; sys.stdout.flush()

# min/max                                                                                             
print 'min(T):', np.min(Ttime37)
print 'max(T):', np.max(Ttime37)

NCL
IDV

IDV is an OPeNDAP tool that can access and display the nature run data. In our OPenDAP server, all files are time aggregated, so they appear as a single dataset for each location.

This is an example to open and display the field 'T' (air temperature) from the collection 'inst01hr_3d_T_Cv'. The OPenDAP URL for this dataset is http://opendap.nccs.nasa.gov:80/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv. The following steps are valid for IDV version 5.0u1 running on a Linux desktop.

From the 'Dashboard' panel

  • Select Data Choosers -> URLS.In the URL field, enter the above OPenDAP URL and click on 'Add Source'
  • Select Field Selector and choose the 3D field'air_temperature'. The 'Times' tab lists all the available levels and times for this data. At this point, one can select specific times, level and regions (subsetting) from the 'Times' and 'Level' and 'Region' tabs. Click on 'Create Display'.

Proprietary clients

Matlab
IDL