Recipe: Matlab program as OPeNDAP client: Difference between revisions

From GEOS-5
Jump to navigation Jump to search
Pchakrab (talk | contribs)
1,048 edits
Created page with "Back to G5NR Data Access Guide. == Problem == By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server http://opendap.nccs.nasa.gov/dods/OSSE/G5..."
 
Move from http to https
 
(5 intermediate revisions by one other user not shown)
Line 4: Line 4:


By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server
By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server
  http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km
  https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km


we want to read the surface temperature data inside the box bound by latitudes 25<sup>o</sup>N, 50<sup>o</sup>N and longitudes -130<sup>o</sup>W, -65<sup>o</sup>W for 2006/Sep/18, 9z and compute its min/max.
we want to read the surface temperature data inside the box bound by latitudes 25<sup>o</sup>N, 50<sup>o</sup>N and longitudes -130<sup>o</sup>W, -65<sup>o</sup>W for 2006/Sep/18, 9z and compute its min/max.


== Solution ==
== Solution ==
First, we ensure that our NetCDF-4 library has been built with Fortran and OPeNDAP support. If both the queries
MATLAB version 2012a and later has native openDAP support. The following code has been tested with MATLAB version 2014a.
> nc-config --has-f90
> nc-config --has-dap
 
return "yes", we have a compatible NetCDF-4 library. Here <code>nc-config</code> is a utility bundled with NetCDF-4 package.


The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at
The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at
  http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info
  https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info


==== Code ====
==== Code ====
This code accesses the collection <code>inst01hr_3d_T_Cv</code> from the OPeNDAP server and reads a subset of the temperature data.
This code accesses the collection <code>inst01hr_3d_T_Cv</code> from the OPeNDAP server and reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min. It then computes the max/min of the above data at the surface (level=72).


NOTE:
NOTE:
# Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
# Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
# While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase. This is reflected in the call to <code>nf90_inq_varid</code>.
# While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase.
# Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly <code>inst</code> files are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772<sup>th</sup> step.
# Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly <code>inst</code> files are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772<sup>th</sup> step.


<syntaxhighlight lang="fortran" line>
<syntaxhighlight lang="matlab" line>
program g5nr_reader_dap
% opendap url
 
opdurl = 'https://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv'
  use netcdf          ! for reading the NR files
 
  implicit none
 
  ! File name
  ! ---------
  character(len=256) :: T_file


  ! 4D array: (lon,lat,lev,time)
% bounding box
   ! ------------------------------------
%   lons = -130:0.5:-65
   real, allocatable :: T(:,:,:,:)
%   lats =  25:0.5:50
imin = round((-130. + 180.)/0.5);
imax = round(( -65. + 180.)/0.5);
jmin = round(( 25 +  90.)/0.5);
jmax = round((  50 +  90.)/0.5);
% corresponding array sizes
im = imax-imin+1;
jm = jmax-jmin+1;
lm = 72; % read all 72 levels


  ! Miscellaneous
% read temperature inside the bounding box
  ! -------------
% time 11772 corresponds to 2006/09/18, 9z
  integer :: ierr
start = [imin, jmin, 1, 11772];
  integer :: im, jm, lm
count = [im, jm, lm, 1];
  integer :: ncid, varid
fprintf('Reading Tsub (subset of T)...');
  integer :: start(4), count(4)
Tsub = ncread(opdurl, 't', start, count);
  integer :: imin, imax, jmin, jmax
fprintf('done.\n')
  real :: minlat, minlon, maxlat, maxlon


  ! We replace filename by opendal url
% compute max/min of Tsub
  ! ----------------------------------
fprintf('max(Tsub): %f\n', max(Tsub(:)));
  T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"
fprintf('min(Tsub): %f\n', min(Tsub(:)));


  ! file name and (subset) dimensions
% compute max/min of Tsub at surface
  ! bounding box:
level = 72; % surface
  !  lons = -130:0.5:-65
Tsub_surface = Tsub(:,:,72);
  !  lats =  25:0.5:50
fprintf('max(Tsub at surface): %f\n', max(Tsub_surface(:)));
  ! ---------------------------------
fprintf('min(Tsub at surface): %f\n', min(Tsub_surface(:)));
  ! array indices from the
  ! bounding lat/lon values
  imin = nint((-130. + 180.)/0.5)
  imax = nint(( -65. + 180.)/0.5)
  jmin = nint((  25 +  90.)/0.5)
  jmax = nint((  50 +  90.)/0.5)
  ! array sizes
  im = imax-imin+1
  jm = jmax-jmin+1
  lm = 72
 
  ! Open url and get var id
  ! -----------------------
  call check(nf90_open(T_file,NF90_NOWRITE,ncid), "opening T file")
  call check(nf90_inq_varid(ncid,"t",varid), "getting T varid")
 
  ! Read temperature data
  ! ---------------------
  allocate(T(im,jm,lm,1))            ! global 4D array with 1 time level
  start = [imin, jmin, 1, 11772]    ! time level 11772 corresponds to 2006/09/18, 9z
  count = [im, jm, lm, 1]            ! 1 time level, 3D (lon,lat,lev) array
  write(*,*) 'Reading T...'
  call check(nf90_get_var(ncid,varid,T,start=start,count=count), "reading T")
  write(*,*)'T: ', maxval(T),minval(T)
 
  ! close file, release memory
  call check(nf90_close(ncid), "closing T file")
  deallocate(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_dap
</syntaxhighlight>
</syntaxhighlight>


==== Compile and link ====
==== Output ====
 
Running this MATLAB script, we get the output
We use the utility <code>nf-config</code> (included in the NetCDF-4 installation) to identify the linking rules
 
For a typical NetCDF-4 installation, the above code, <code>g5nr_reader_dap.f90</code> can be compiled and linked to the NetCDF-4 library via
gfortran -o g5nr_reader_dap.x `nf-config --fflags` g5nr_reader_dap.f90 `nf-config --flibs`
 
creating the executable <code>g5nr_reader_dap.x</code>.
 
NOTE:
# You can use your favorite Fortran compiler instead of <code>gfortran</code>.
# If the NetCDF-4 library was built with parallel I/O support, you will need to use <code>mpif90</code> to link, even if your code does not use the MPI library.
 
==== Run ====
Running the executable,
 
./g5nr_reader_dap.x


produces the output
>> g5nr_reader
  Reading T (subset)...
  T (subset):   305.6048      191.6956
opdurl =
https://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv
  Reading Tsub (subset of T)...done.
  max(Tsub): 305.604828
min(Tsub): 191.695648
max(Tsub at surface): 305.604828
min(Tsub at surface): 276.292328


== Discussion ==
== Discussion ==

Latest revision as of 10:25, 9 April 2019

Back to G5NR Data Access Guide.

Problem

By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server 

https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km

we want to read the surface temperature data inside the box bound by latitudes 25oN, 50oN and longitudes -130oW, -65oW for 2006/Sep/18, 9z and compute its min/max.

Solution

MATLAB version 2012a and later has native openDAP support. The following code has been tested with MATLAB version 2014a.

The metadata for the collection inst01hr_3d_T_Cv is available at 

https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info

Code

This code accesses the collection inst01hr_3d_T_Cv from the OPeNDAP server and reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min. It then computes the max/min of the above data at the surface (level=72).

NOTE:

  1. Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
  2. While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase.
  3. Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly inst files are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772th step.
% opendap url
opdurl = 'https://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv'

% bounding box
%   lons = -130:0.5:-65
%   lats =   25:0.5:50
imin = round((-130. + 180.)/0.5);
imax = round(( -65. + 180.)/0.5);
jmin = round((  25 +  90.)/0.5);
jmax = round((  50 +  90.)/0.5);
% corresponding array sizes
im = imax-imin+1;
jm = jmax-jmin+1;
lm = 72; % read all 72 levels

% read temperature inside the bounding box
% time 11772 corresponds to 2006/09/18, 9z
start = [imin, jmin, 1, 11772];
count = [im, jm, lm, 1];
fprintf('Reading Tsub (subset of T)...');
Tsub = ncread(opdurl, 't', start, count);
fprintf('done.\n')

% compute max/min of Tsub
fprintf('max(Tsub): %f\n', max(Tsub(:)));
fprintf('min(Tsub): %f\n', min(Tsub(:)));

% compute max/min of Tsub at surface
level = 72; % surface
Tsub_surface = Tsub(:,:,72);
fprintf('max(Tsub at surface): %f\n', max(Tsub_surface(:)));
fprintf('min(Tsub at surface): %f\n', min(Tsub_surface(:)));

Output

Running this MATLAB script, we get the output 

>> g5nr_reader

opdurl =

https://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv

Reading Tsub (subset of T)...done.
max(Tsub): 305.604828
min(Tsub): 191.695648
max(Tsub at surface): 305.604828
min(Tsub at surface): 276.292328

Discussion

See Also

  1. File Spec: File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf
  2. Recipe: Fortran program to read data from downloaded file

No Warranty

Copyright

Retrieved from "https://geos5.org/wiki/index.php?title=Recipe:_Matlab_program_as_OPeNDAP_client&oldid=4304"