How to define geographicGrid in matlab code?

I am using following matlab code to extract data but it is giving some error.

img_file = 'c:\data\4MARCH24_BAND7_SUBSET.dat'
img = hypercube(img_file);
shapefile = 'c:\data\Achi Khurd.shp';
S = shaperead(shapefile);

%remove small triangles and squares
coord_lengths = arrayfun(@(s) numel(s.X), S);
not_useful_mask = coord_lengths < 10;
S(not_useful_mask) = [];
polygon = polyshape({S.X}, {S.Y});

% Create a logical mask
[Z,R] = readgeoraster(img_file);
[lat,lon] = geographicGrid(R);
[Lon, Lat] = meshgrid(lon, lat);
logical_mask = reshape(isinterior(polygon, Lon(:), Lat(:)), size(Lon)); 

% Use the logical mask to extract data 
extracted_data = img(logical_mask);
disp(extracted_data);

It is giving following error

Incorrect number or types of inputs or outputs for function geographicGrid.

Error in Test (line 14)

[lat,lon] = geographicGrid(R);

I am attaching the data also in order to run the above mentioned matlab code.

I would be thankful for suggesting me to fix it.

Deva

Answer 1

[Cris LaPierre]

Ultimately, the issue is because your input to geographicGrid is the wrong data type.

geographicGrid expects the input to be a GeographicCellsReference or GeographicPostingsReference object (see here). However, if you inspect the class of your variable R in the Workspace, you will see it is a MapCellsReference object.

Basically, your data file does not contain the data format you expect.

You can compare the difference between the two as well as see their corresponding object functions on their respective doc pages.

• GeographicCellsReference
• MapCellsReference

[Devendra]

Thanks for your help.

geographicGrid expects the input to be a GeographicCellsReference or GeographicPostingsReference object (see here).

However, if you inspect the class of your variable R in the Workspace, you will see it is a MapCellsReference object.

How do I change format of my data which is MapCellsReference object into GeographicCellsReference or GeographicPostingsReference object ?

Deva

[Cris LaPierre]

This is not an area I am an expert in so could not provide an authoritative answer. How is your data obtained? Do you have a dataset you can share where your code does work?

[Cris LaPierre]

Your shapefile is also in map coordinates, not lat/lon. Given that, I'm not sure you want to convert. Why not use functions that do work witih mapCellsReference objects?

• worldGrid
• mapshow

[Devendra]

I have downloaded the sentinal2 data and subset the data as per shape file using ENVI software. I am attaching the data.

I would appreciate any further suggestions/help in reformating my data into GeographicCellsReference or GeographicPostingsReference object.

Deva

[Devendra]

Thanks for your suggestion. I have used the worldGrid as follows

[Z,R] = readgeoraster(img_file);

p = R.ProjectedCRS;

[x,y] = worldGrid(R);

[lat,lon] = projinv(p,x,y);

%[lat,lon] = geographicGrid(R);

[Lon, Lat] = meshgrid(lon, lat);

However, some different errors have occurred as follows

Error using repmat

Requested 3596966x3596966 (96396.8GB) array exceeds maximum array size preference (15.8GB). This might cause

MATLAB to become unresponsive.

Error in meshgrid (line 61)

xx = repmat(xrow,size(ycol));

Error in Test (line 18)

[Lon, Lat] = meshgrid(lon, lat);

Related documentation

Please suggest me how to move further.

Deva

[Cris LaPierre]

Have you inspected your variables? You likely do not need to use meshgrid. Your lon and lat variables are already matrices.

[Devendra]

Thank you very much for your kind help. It worked. However, it is not displaying the extracted data because of following part of code;

S = shaperead(shapefile);

%remove small triangles and squares

coord_lengths = arrayfun(@(s) numel(s.X), S);

not_useful_mask = coord_lengths < 10;

S(not_useful_mask) = [];

polygon = polyshape({S.X}, {S.Y});

Since it is showing the size of polygon 1 1

in place of actual size of shape file which is showing as 1 22.

I want to see the extracted data from main file for the shape file.

Please suggest me how to get extracted data using shape file.

Deva

[Cris LaPierre]

For the data you have shared, none of the values from Z are within the shapefile boundaries, so the result is correctly a 1x1 empty polygon. Use mapshow to visualize your data.

[Devendra]

I am very sorry for my mistake not to overlay the shape file boundary over image. However, I am again sending the data which contains the shape file boundary over image. I am also sending PNG file showing the overlay of shape file over image. However, matlab code has not extrected the data as per shape file.

I request you to please have a look on this new data set and suggest me how to get the data using shape file.

I really appreciate your kind help.

Deva

[Cris LaPierre]

Your shapefile and dat file are both in map coordinates. You need to dig into your files to see why they do not overlap.

I came to the conclusion that this is because your dat file and shapefile are using different projections.

% Shapefile projections

info = shapeinfo(shapefile);

crs = info.CoordinateReferenceSystem

crs =

projcrs with properties:

Name: "LCC_WGS"

GeographicCRS: [1×1 geocrs]

ProjectionMethod: "Lambert Conic Conformal (2SP)"

LengthUnit: "meter"

ProjectionParameters: [1×1 map.crs.ProjectionParameters]

% dat file projections

R.ProjectedCRS

ans =

projcrs with properties:

Name: "WGS 84 / UTM zone 43N"

GeographicCRS: [1×1 geocrs]

ProjectionMethod: "Transverse Mercator"

LengthUnit: "meter"

ProjectionParameters: [1×1 map.crs.ProjectionParameters]

You therefore need to use projinv on both your files.

S = shaperead(shapefile);
[S.lat,S.lon] = projinv(crs,S.X,S.Y);


[Z,R,cmap] = readgeoraster(img_file);
Z = double(Z);
p = R.ProjectedCRS;
[x,y] = worldGrid(R);
[lat,lon] = projinv(p,x,y);

Then at last, the 2 shapes overlap, and your mask should no longer be empty.

[Devendra]

Thanks a lot for your valuable suggestions. However, it is giving following error;

Unrecognized function or variable 'crs'.

Error in Test (line 9)

[S.lon,S.lat] = projinv(crs,S.X,S.Y);

Deva

[Voss]

crs is defined in Cris's code above:

% Shapefile projections
info = shapeinfo(shapefile);
crs = info.CoordinateReferenceSystem;

[Devendra]

Thank you very much. It worked and extracted data from single band image when I use

polygon = polyshape({S.lon}, {S.lat});

in place of

polygon = polyshape({S.X}, {S.Y});

However, when I wanted to extract the data from 7 band image it does not extract any data. I request you to please help me out in fixing out this error also.

Deva

[Cris LaPierre]

I don't have enough informatino to offer help. Again, this is not a space I'm an expert in. I'm just reading doc pages. I would start by looking at the variable containing a single band image and compare it to the one of a 7 band image. How is the variable different? Does the code you've written work the same on both variiables? If not, what changes to the code (or variable) can you make to ensure it works for both?

[Devendra]

Please ignore my above mentioned request. It has already been done. I am very grateful to you for your help.

Deva

[Devendra]

I am using two different shape files of same type. However, one shape file abdualpur is not giving any error while second shape file is giving following error

BAGHPAT_VILLAGES_RECT.shp

Error using projcrs/projinv

Too many input arguments.

Error in ds_extract_data (line 21)

[S.lon,S.lat] = projinv(crs,S.X,S.Y);

I am attaching both shape files and request you to please have a look on them and suggest me how to use BAGHPAT_VILLAGES_RECT.shp without error.

Devendra

[Cris LaPierre]

BAGHPAT_VILLAGES_RECT.shp contains two polygons. You have to tell projinv which one to use. Use intexing to extract one polygon. Here's an example.

[S(1).lon,S(1).lat] = projinv(crs,S(1).X,S(1).Y);

[Devendra]

Thanks a lot for your help, It worked well.

Devendra

[Devendra]

I am bit hesitent to request you little more help. Actually I have very recently switched on to MATLAB and therefore, I am in a learning phase. My humble request to you to please suggest me how to use customize enviwrite matlab code (attached) for the map_info (attached) which is generated by satellite data in jp2 format to make envi header file similar to envi headre file (attached) . I am attaching these three files and request your kind cooperation.

Devendra

[Cris LaPierre]

I would suggest starting a new question for this. That will give the broader community an opportunity to answer as well, rather than just me.

[Devendra]

Thank you so much for your kind cooperation and suggestion. Certainly I will put forth this question as a new question to larger matlab community. Once again thank you very much for the help. Devendra

[Devendra]

I have put across the question of customizing enviwrite code to larger MATLAB community but so far nobody has responded in last two days. Actually I want to customize enviwrite matlab code for reading my map_info and produce ENVI header file containing map_info and projection_info. Please suggest me whom should I approach to get it coustomized for my purpose. I would appreciate your kind suggestions and would be highly obliged to you.

Devendra

[Cris LaPierre]

Sorry. It can take a while, especially if the topic is niche. I'm afraid this is not an area where I consider myself an expert. I do not know who to direct you to, and I am unfamiliar with the topics you are asking about.