Spatial subset of geotiff through masking by shapefile?

40 views (last 30 days)
I would like to create a mask (using poly2mask) to apply to a geotiff image, based on a single polygon in shapefile .shp format.
I have put together simple code below but evidently I have something wrong with syntax (maybe even just totally improper usage) of the worldtoIntrinsic part.
I have searched help avenues, but haven't cracked it.
Error message: (An issue with the structure of R, that I am returning by geotiffread?)
Undefined function 'worldToIntrinsic' for input arguments of type 'map.rasterref.GeographicPostingsReference'. Error in polygon_shapefileread_example (line 15) [ix, iy] = worldToIntrinsic(R,rx,ry);
Attached are the small geotiff and shapefile.
Many thanks for any help, Matt
%Read in GeoTIFF
[I R] = geotiffread('geotiff_example.tif');
% Read Region of Interest shapefile
roi = shaperead('shapefile_example.shp');
% Remove trailing nan from shapefile
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);
% convert to image coordinates
[ix, iy] = worldToIntrinsic(R,rx,ry);
%Make the mask
mask = poly2mask(ix,iy,R.RasterSize(1),R.RasterSize(2));
%Following line checks some 1's are generated in mask
maskcheck=sum(sum(mask));

Accepted Answer

Chad Greene
Chad Greene on 4 Mar 2015
Does this work for you?
% Get image info:
R = geotiffinfo('geotiff_example.tif');
% Get x,y locations of pixels:
[x,y] = pixcenters(R);
% Convert x,y arrays to grid:
[X,Y] = meshgrid(x,y);
roi = shaperead('shapefile_example.shp');
% Remove trailing nan from shapefile
rx = roi.X(1:end-1);
ry = roi.Y(1:end-1);
mask = inpolygon(X,Y,rx,ry);
  11 Comments
kinga kulesza
kinga kulesza on 10 Nov 2022
for multiple polygons:
[X,Y]=meshgrid(lon,lat);
roi = shaperead('multiple_polygons.shp');
for i=1:numel(roi)
rx = roi(i).X(1:end-1);
ry = roi(i).Y(1:end-1);
mask(:,:,i) = inpolygon(X,Y,rx,ry);
end

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!