How to change the given code for UTM projection?

I am using the following code to extract the data using shape file. % Load the georeferenced raster and its spatial reference object [img, R] = readgeoraster(img_file); % Load the shapefile and retrieve its info S = shaperead(shape_file); info = shapeinfo(shape_file); % Retrieve the Coordinate Reference System from the shapefile info crs = info.CoordinateReferenceSystem; % Transform the coordinates of the second shape in the shapefile from projected to geographic [S(2).lat, S(2).lon] = projinv(crs, S(2).X, S(2).Y); % Open a figure window and set the basemap to satellite imagery figure; geobasemap('satellite'); % Display the second shape on the satellite basemap geoplot(S(2).lat, S(2).lon); % Generate a geographic raster grid projectedCRS = R.ProjectedCRS; [x, y] = worldGrid®; [lat, lon] = projinv(projectedCRS, x, y); % Construct a polygon from the second shape's geographic coordinates polygon = polyshape({S(2).lon}, {S(2).lat}); % Generate a mask by testing if the raster's geographic grid points are inside the polygon mask = reshape(isinterior(polygon, lon(:), lat(:)), [length(lat), length(lon)]);

Actually I am interested in UTM projection instead of geographic projection as done in above mentioned code. Since the source image jp2 is in UTM projection and I want to keep the same projection in extracted image also. May I request you to please suggest me how to use it for UTM projection in place of geographic projection?

Answer 1

[Vinayak]

Hi Gauri,

A few suggestions for your next query: prefer attaching the code as code text to make it more readable. If there are multiple files or a large function, you can also attach files using the "Paper clip" icon in the insert section.

Regarding your code, you only need the geo-coordinates for use with geoplot. You can save them as extra variables instead of replacing the default UTM coordinates. This way, you can generate the worldGrid using the default UTM coordinates.

Here's an example:

[lat, lon] = projinv(crs, S(2).X, S(2).Y);
geoplot(lat, lon);
[x, y] = worldGrid(R);

This approach should eliminate the need for unnecessary transformations with the projected CRS.