You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to define geographicGrid in matlab code?
11 views (last 30 days)
Show older comments
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
Accepted Answer
Cris LaPierre
on 29 Mar 2024
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.
26 Comments
Devendra
on 29 Mar 2024
Edited: Devendra
on 29 Mar 2024
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
on 29 Mar 2024
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
on 29 Mar 2024
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
on 29 Mar 2024
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
on 30 Mar 2024
Have you inspected your variables? You likely do not need to use meshgrid. Your lon and lat variables are already matrices.
Devendra
on 30 Mar 2024
Edited: Devendra
on 30 Mar 2024
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
on 30 Mar 2024
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
on 30 Mar 2024
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
on 30 Mar 2024
Edited: Cris LaPierre
on 30 Mar 2024
Your shapefile and dat file are both in map coordinates. You need to dig into your files to see why they do not overlap.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1655901/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1655906/image.png)
Devendra
on 30 Mar 2024
Edited: Devendra
on 30 Mar 2024
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
on 31 Mar 2024
Edited: Cris LaPierre
on 1 Apr 2024
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
on 20 Apr 2024
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
on 21 Apr 2024
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
on 22 Apr 2024
Edited: Devendra
on 22 Apr 2024
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
on 22 Apr 2024
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
on 24 Apr 2024
Edited: Devendra
on 24 Apr 2024
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
on 24 Apr 2024
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.
Aksh
on 27 May 2024
>> indices_caf
2481 2259
Ale.shp
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.
> In polyshape/checkAndSimplify (line 526)
In polyshape (line 175)
In ds_hants_indices_caf (line 45)
Requested 408x5604579 (17.0GB) array exceeds maximum array size preference (15.8GB). This might cause MATLAB to become unresponsive.
Error in polyshape/isinterior>check_inpolygon (line 65)
x = x(ones(Nv,1),:);
Error in polyshape/isinterior (line 50)
[in, on] = check_inpolygon(X, Y, xv, yv, tol);
Error in indices_caf (line 53)
mask = reshape(isinterior(polygon, lon(:), lat(:)), size(lon));
May I request you to please suggest me how to resolve this error.
I will be thankful to you.
Cris LaPierre
on 27 May 2024
Please ask a new question as this one is already answered. When doing so, please provide enough of the code for us to reporduce the error. It would be helpful to attach the necessary file(s) to your post using the paperclip icon. You may have to zip them first.
More Answers (0)
See Also
Categories
Find more on Data Type Identification in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)