Issue: H5 dataset extraction based on shapefile and Datetime

7 visualizaciones (últimos 30 días)
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 comentarios
Anmol Dhiman
Anmol Dhiman el 27 de En. de 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 el 27 de En. de 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,[]);

Iniciar sesión para comentar.

Respuesta aceptada

Anmol Dhiman
Anmol Dhiman el 31 de En. de 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 comentario
Muhammad Usman Liaqat
Muhammad Usman Liaqat el 3 de Feb. de 2021
Dear Anmol,
Thanks for your reply but I found alternative solution to resolve it.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Preprocessing Data en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by