Prog_Name = 'ReadGrid.m';
FILE_NAME = 'AMSR_E_L3_RainGrid_B05_200707.hdf';
GRID_NAME = 'MonthlyRainTotal_GeoGrid'
DATAFIELD_NAME = 'RrLandRain'

%Open the HDF-EOS2 Grid File
file_id = hdfgd('open', FILE_NAME, 'rdonly')
if (file_id == -1)
   disp([Prog_Name ': ERROR - ' FILE_NAME ' could not be opened.'])
   return; end

%Open a grid named 'MonthlyRainTotal_GeoGrid'
grid_id = hdfgd('attach', file_id, GRID_NAME)
if (grid_id == -1)
   disp([Prog_Name ': ERROR - Unable to open grid ' GRID_NAME])
   status = hdfgd('close', file_id);
   return; end

%Get information about the spatial extents of the grid.
[xdimsize, ydimsize, upleft, lowright, status] = hdfgd('gridinfo', grid_id)
if (status == -1)
   hdfgd('close', file_id);
   disp([Prog_Name ': ERROR - Unable to open grid' GRID_NAME])
   return; end

%Reading Data from a Data Field
[rrland, fail] = hdfgd('readfield', grid_id, DATAFIELD_NAME, [], [], []);
if (fail == -1)
   disp([Prog_Name ': ERROR - Unable to read data from' DATAFIELD_NAME])
   status = hdfgd('detach', grid_id);
   status = hdfgd('close', file_id);
   return; end

%We need to readjust the limits of latitude and longitude. HDF-EOS is using DMS(DDDMMMSSS.SS) format to represent degrees.
%So to calculate the lat and lon in degree, one needs to convert minutes and seconds into degrees. 

%The following is the detailed description on how to calculate the latitude and longitude range based on lowright and upleft.
% One should observe the fact that 1 minute is 60 seconds and 1 degree is 60 minutes. 

%First calculate the difference of .SS between lowright and upleft:
lowright_ss= lowright*100-floor(lowright)*100;
upleft_ss = upleft*100-floor(upleft)*100;
dss = lowright_ss - upleft_ss;

%Then calculate the difference of SSS between lowright and upleft:
lowright_s = mod(floor(lowright),1000);
upleft_s = mod(floor(upleft),1000);

ds =lowright_s - upleft_s +dss/100;

%Then calculate the difference of MMM between lowright and upleft:

lowright_m = mod(floor(lowright/1000),1000);
upleft_m = mod(floor(upleft/1000),1000);

dm = lowright_m-upleft_m +ds/60;

%Then calculate the difference of DDD between lowright and upleft:
lowright_d = floor(lowright/1000000);
upleft_d = floor(upleft/1000000);
dd = lowright_d-upleft_d+dm/60;


lat_limit = dd(2);
lon_limit = dd(1);


%We need to calculate the grid space interval between two adjacent points

scaleX = lon_limit/xdimsize;
scaleY = lat_limit/ydimsize;

%By default HDFE_CENTER is assumed for the offset value, which assigns 0.5 to both offsetX and offsetY.
offsetX = 0.5;
offsetY = 0.5;

%Since this grid is using geographic projection, the latitude and longitude value will be calculated based on the formula:
% (i+offsetX)*scaleX+leftX  for longitude and (i+offsetY)*scaleY+leftY for latitude.


for i = 0:(xdimsize-1)
  lon_value(i+1) = (i+offsetX)*(scaleX) + upleft_d(1);
end

for j = 0:(ydimsize-1)
  lat_value(j+1) = (j+offsetY)*(scaleY) + upleft_d(2);
end

% To use Matlab's geoshow function, in this example,
% the longitude needs to be translated from 0 to 360 degree to -180 to 180 degree,
% the latitude needs to be translated from "north to south"( latitude decreases) to  "south to north"(latitude increases).

lon_offset = -180;
lon_value = lon_offset + lon_value;
lat_value = fliplr(lat_value);

% The following shows the size and order of lon,lat and rrland so that one can see the need of rrland to be transposed.
% This is kind of "help information" for the next operation.
lon_size=size(lon_value);
lat_size=size(lat_value);
rrland_size=size(rrland);
disp('lon size');
disp(lon_size);
disp('lat_size');
disp(lat_size);

%To ensure that the contour overlays with the world map correctly, the fast changing dimension of the plotted data field(rrland) 
%must be the same size as the geo-location dimension specified by the first argument of the contour function. 
%In our case, the first argument represents longitude. the data field needs to be transposed.
ts=transpose(rrland);

%Since we want to plot the data from -180 to 180 in longitude, and the original data starts from 0 to 360(really means 0 to 180, -180 to 0 in longitude),
%So we need to move the data from -180 to 0 before 0 to 180. That means switching the second half of the array with the first half of the array along
% the dimension of longitude.

halfx = floor(xdimsize/2);

ts_reverse=[ts(:, (halfx+1):xdimsize) ts(:, 1:halfx)];

%MATLAB wants to plot the data using lower left as origin of the data. 
%However, as we see from the information obtained by 'gridinfo', the origin is upper left for the file. Hence we need to flip the data.
data = flipud(ts_reverse);

%Visualizing the data field on a world map
contour(lon_value, lat_value, data)
geoshow('landareas.shp', 'FaceColor', [0.4 0.4 0.4])
colorbar

%Dumping a subset of the values in the data field
startvector = [0 0]
stride = []
endvector = [5 5]
[rrland_subset, subset_fail] = hdfgd('readfield', grid_id, DATAFIELD_NAME, startvector, stride, endvector);
if (subset_fail == -1)
   disp([Prog_Name ': ERROR - Unable to read and dump a subset of data from' DATAFIELD_NAME])
   status = hdfgd('detach', grid_id);
   status = hdfgd('close', file_id);
   return; end

%Alternative method to dump the subset of the values in the data field
for i = 1:5
   for j = 1:7
      rrland(i,j);
   end
end

%Close the grid named 'MonthlyRainTotal_GeoGrid'
hdfgd('detach', grid_id);

%Close the HDF-EOS2 Grid File
hdfgd('close', file_id);