How to define geographicGrid in matlab code?

10 次查看（过去 30 天）

显示更早的评论

Devendra 2024-3-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2100491-how-to-define-geographicgrid-in-matlab-code

评论： Cris LaPierre 2024-5-27

采纳的回答： Cris LaPierre

Data.zip

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Cris LaPierre 2024-3-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2100491-how-to-define-geographicgrid-in-matlab-code#answer_1433421

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 个评论
显示 24更早的评论隐藏 24更早的评论

Cris LaPierre 2024-3-30

编辑：Cris LaPierre 2024-3-30

在 MATLAB Online 中打开

Data.zip

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.

Aksh 2024-5-27

>> 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 2024-5-27

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.

请先登录，再进行评论。