이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How to define geographicGrid in matlab code?
조회 수: 6 (최근 30일)
이전 댓글 표시
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
채택된 답변
Cris LaPierre
2024년 3월 29일
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
Devendra
2024년 3월 29일
편집: Devendra
2024년 3월 29일
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
2024년 3월 29일
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
2024년 3월 29일
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
2024년 3월 29일
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
2024년 3월 29일
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
2024년 3월 30일
Have you inspected your variables? You likely do not need to use meshgrid. Your lon and lat variables are already matrices.
Devendra
2024년 3월 30일
편집: Devendra
2024년 3월 30일
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
2024년 3월 30일
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
2024년 3월 30일
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
2024년 3월 30일
편집: Cris LaPierre
2024년 3월 30일
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.
Voss
2024년 3월 30일
crs is defined in Cris's code above:
% Shapefile projections info = shapeinfo(shapefile); crs = info.CoordinateReferenceSystem;
Devendra
2024년 3월 30일
편집: Devendra
2024년 3월 30일
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
2024년 3월 31일
편집: Cris LaPierre
2024년 4월 1일
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
2024년 4월 1일
Please ignore my above mentioned request. It has already been done. I am very grateful to you for your help.
Deva
Devendra
2024년 4월 20일
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
2024년 4월 21일
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
2024년 4월 22일
편집: Devendra
2024년 4월 22일
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
2024년 4월 22일
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
2024년 4월 22일
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
2024년 4월 24일
편집: Devendra
2024년 4월 24일
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
2024년 4월 24일
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.
gauri
2024년 5월 4일
Hello Cris
I have asked a question and link is as follows;
May I request you to please have a look on it and suggest me how to do it?
Gauri
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.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Data Type Identification에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- 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)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)