Issue: H5 dataset extraction based on shapefile and Datetime

6 vues (au cours des 30 derniers jours)
Muhammad Usman Liaqat
Muhammad Usman Liaqat le 23 Jan 2021
I am trying to seprate data from H5 file, reporject and subset it based on shapefile. My data is available in DateTime Format. However, I am facing error while trying to seprate based on DateTime.
%load shapefiles of indus subbasin lat/lon boundaires
indusFile='UIB_Subbasubs.shp';
S=shaperead(indusFile);
info=shapeinfo(indusFile);
%X Y data is Lon Lat data - shapefile is geographic, not projected
[S.Lon]=S.X;
[S.Lat]=S.Y;
S = rmfield(S,{'X','Y'});
% Indus subbasins
indusBasins = {'Astore';'Gilgit';'Hunza';'Indus Downstream';'Shyok';'Shingo';'Shigar';' Zanskar'};
% union of subbasins to clip parbal output to
%(could extract by individual subbasins instead)
[idxSN,~] = ismember({S.Name},indusBasins);
tempIdx=find(idxSN);
for i = 1:length(indusBasins)
geoin = polyshape(S(tempIdx(i)).Lon,S(tempIdx(i)).Lat);
% combine masks from each subbasin
if i==1
basinUnion=geoin;
else
basinUnion=union(basinUnion,geoin);
end
end
%yr = 2000;
sweDateStart = datetime(2000,01,01,0,0,0);
sweDateEnd = datetime(2000,01,01,23,0,0);
sweDates = sweDateStart:hours(1):sweDateEnd;
sweDates.Format = 'yyyyMMddHHmmSS';
%example parbal output from one year
sweFile= ['D:/matlab tutorial/MODIS/20001001_Indus463m_sin_forcings.h5'];
for d = 1:length(sweDates)
% sweDate=datenum(2006,1,1);
[x,hdr,h5mdates]=getMelt(sweFile,'presZ',sweDates(d));
%get geographic raster reference object for the sinusodial data grid
[ geoSWE, geoSWE_RefMatrix, geoSWE_RasterReference] = ...
rasterReprojection(x,hdr.RasterReference,hdr.ProjectionStructure,[]);%,varargin )
% extract SWE for each subbasin
[indusBoundary, indusR] = vec2mtx(basinUnion.Vertices(:,2),...
basinUnion.Vertices(:,1), geoSWE, geoSWE_RefMatrix,'filled');
assert(isequal(geoSWE_RefMatrix,indusR),'refmat issues...')
%0==inside basin boundary, 1== intersects basin boundary, 2== outside basin
%boundary
validPxls=indusBoundary~=2;
basinSWE=geoSWE;
basinSWE(~validPxls)=nan;
% %save geographic referenced SWE output
% saveFile='F:\matlab tutorial\MODIS\',d ,'.tif';
% geotiffwrite(saveFile,basinSWE,geoSWE_RefMatrix);
D=num2str(d)
saveFile= strcat(D,'SWE_DOY.tif')
geotiffwrite(saveFile,basinSWE,geoSWE_RefMatrix);
%check if it worked
info=geotiffinfo(saveFile);
end
  2 commentaires
Anmol Dhiman
Anmol Dhiman le 27 Jan 2021
Hi Usman,
Could you share the part of input data as well for better debugging of the issue.
Muhammad Usman Liaqat
Muhammad Usman Liaqat le 27 Jan 2021
Dear Anmol,
Kindly find Input data from following link. My overall aim to project h5 from sinusoidal to wgs84 projection and save in GeoTiff format.
%I also tried to use this data in simple code by applying following function
fname='20001001_Indus463m_sin_forcings.h5';
hdr=GetCoordinateInfo(fname,'/Grid',[1702 3636]);
presZ=h5read(fname,'/Grid/presZ');
[presZ_geog,R,RR]=rasterReprojection(presZ,hdr.RefMatrix,hdr.ProjectionStructure,[]);

Connectez-vous pour commenter.

Réponse acceptée

Anmol Dhiman
Anmol Dhiman le 31 Jan 2021
Hi Usman,
The above error is coming because the formats you are comparing are not equal.
In line number 36, getMelt
h5mdates =
24×1 string array
and
matdates =
datetime
To avoid above error, use
find(datestr(matdates,'yyyyMMddHHmmSS')==h5mdates)
You also need to add a condition in case the particular data is not found before executing line 37.
Hope it Helps
  1 commentaire
Muhammad Usman Liaqat
Muhammad Usman Liaqat le 3 Fév 2021
Dear Anmol,
Thanks for your reply but I found alternative solution to resolve it.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Preprocessing Data dans Help Center et File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by