GEOS-5

Recipe: Matlab program as OPeNDAP client: Difference between revisions

Newer edit →
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..."
 
Pchakrab (talk | contribs)
1,048 edits
Line 9: Line 9:


== 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
Line 23: Line 19:
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 = 'http://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)
  ! ------------------------------------
  real, allocatable :: T(:,:,:,:)
 
  ! Miscellaneous
  ! -------------
  integer :: ierr
  integer :: im, jm, lm
  integer :: ncid, varid
  integer :: start(4), count(4)
  integer :: imin, imax, jmin, jmax
  real :: minlat, minlon, maxlat, maxlon
 
  ! We replace filename by opendal url
  ! ----------------------------------
  T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"
 
  ! file name and (subset) dimensions
  ! bounding box:
  !  lons = -130:0.5:-65
  !  lats =  25:0.5:50
  ! ---------------------------------
  ! 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
% bounding box
    character(len=*), intent(in) :: loc
%  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


    if(status /= NF90_NOERR) then
% read temperature inside the bounding box
      write (*,*) "Error at ", loc
% time 11772 corresponds to 2006/09/18, 9z
      write (*,*) NF90_STRERROR(status)
start = [imin, jmin, 1, 11772];
    end if
count = [im, jm, lm, 1];
fprintf('Reading Tsub (subset of T)...');
Tsub = ncread(opdurl, 't', start, count);
fprintf('done.\n')


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


end program g5nr_reader_dap
% 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(:)));
</syntaxhighlight>
</syntaxhighlight>


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