Issue: H5 dataset extraction based on shapefile and Datetime

Question

Muhammad Usman Liaqat 2021-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/724897-issue-h5-dataset-extraction-based-on-shapefile-and-datetime

评论： Muhammad Usman Liaqat 2021-2-3

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 2021-1-27

Hi Usman,

Could you share the part of input data as well for better debugging of the issue.

Muhammad Usman Liaqat 2021-1-27

在 MATLAB Online 中打开

code.png

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.

https://drive.google.com/file/d/1BfqE3hwgLLeJcSC2Hj8cxCv5zLXuwlyW/view?usp=sharing

%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,[]);