フィルターのクリア

Issue: H5 dataset extraction based on shapefile and Datetime

3 ビュー (過去 30 日間)
Muhammad Usman Liaqat
Muhammad Usman Liaqat 2021 年 1 月 23 日
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 件のコメント
Anmol Dhiman
Anmol Dhiman 2021 年 1 月 27 日
Hi Usman,
Could you share the part of input data as well for better debugging of the issue.
Muhammad Usman Liaqat
Muhammad Usman Liaqat 2021 年 1 月 27 日
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,[]);

サインインしてコメントする。

採用された回答

Anmol Dhiman
Anmol Dhiman 2021 年 1 月 31 日
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 件のコメント
Muhammad Usman Liaqat
Muhammad Usman Liaqat 2021 年 2 月 3 日
Dear Anmol,
Thanks for your reply but I found alternative solution to resolve it.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangePreprocessing Data についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by