Issue when masking a raster map by a shape file using 'vec2mtx'

Question

0 个投票

Hi, I'm trying to mask a raster map by a shape file. The problem is at the end the masked region is shifted. I think it should be something wrong with 'vec2mtx'. Does anyone have a suggestion on this? The code is simply:

    roi = shaperead(shape_name);
    info = geotiffinfo(map_name);
    [A, R] = geotiffread(map_name);
    zone = '38S';
    [Latlim, Lonlim] = projinv(info, R.XWorldLimits, R.YWorldLimits);
    refmat = makerefmat('RasterSize', R.RasterSize, 'Latlim', Latlim, ...
      'Lonlim', Lonlim, 'ColumnsStartFrom', R.ColumnsStartFrom, ...
      'RowsStartFrom', R.RowsStartFrom);
    rx = roi.X(1:end-1);
    ry = roi.Y(1:end-1);
    [shape_lat, shape_lon] = projinv(info, rx, ry);
    [A1, R2] = vec2mtx(shape_lat, shape_lon, A, refmat, 'filled');
    fin_idx = (A1 ~= 2);
    A2 = A;
    A2(fin_idx) = NaN;
    geotiffwrite(out_map_name, A2, R,...
      'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);