How to create a map mask using coordinates

21 次查看(过去 30 天)
Hello,
I have data (temperature) and their coordinates associated. I have two vectors with latitude and longitude defining an area, let say it's a square for easiness (in fact I have thousands of coordinates) having the coordinates :
latX = [10.7; 10.7; 20.2; 20.2];
lonX = [-15.3; -5.1; -5.1; -15.3];
And X represents some data I have access to for the whole area.
My coordinates have the from :
lat_t = 0:1:30;
lon_t = -20:1:0;
[lon,lat]=meshgrid(lon_t,lat_t);
Then I apply the poly2mask to return a logical value of 1 for points inside the square:
bw = poly2mask(latX,lonX,size(lat,1),size(lat,2));
It doesn't work because I have negative coordinates, so it can not use the indexs. What could the solution be ? Is there a poly2mask method dedicated to coordinates ?

采纳的回答

MaHa
MaHa 2021-8-9
One solution I found for those who will have the same problem. I used knnsearch to find the closest points in my latitude and longitude grid, and used their indexs and applied poly2mask on it.
clear indLat indLon latPix lonPix coord
for i = 1:size(area,1)
indLat = knnsearch(lat(1,:)',area(i,2));
indLon = knnsearch(lon(:,1),area(i,1));
latCoord(i,1) = indLat;
lonCoord(i,1) = indLon;
end
bw = poly2mask(latCoord(:,1),lonCoord(:,1),size(lon,1),size(lon,2));
with area being my matrix of real coordinates associated to the area (column 1 = lat, 2 = lon). I can use this solution because my latitude and longitude grid are of high resolution enough, and that I am not constrained by the quality of the final results, it is ok if I miss a pixel, which may not be for very high resolution tasks.
  1 个评论
Sean de Wolski
Sean de Wolski 2021-8-9
this assumes that latitude and longitude are equal which they're not and so this will give wrong answers.

请先登录,再进行评论。

更多回答(1 个)

Sean de Wolski
Sean de Wolski 2021-8-6
What is your end goal?
You need to project lat and lon into cartesian coordinates. From there, you can use poly2mask. The resulting mask will still be in projected coordinates. You can then inverse-project it back to lat/lon.
Look at projfwd, projinv - you'll need to find an appropriate projection based on your goal (equal area, etc.)
  4 个评论
MaHa
MaHa 2021-8-10
Hum, I don't know because right now it works fine and returns me the answer I was expecting. But you are right I cross 30 lat degrees and my pixels are not of equal area. I have a regular grid (in degrees) but not in distance.
Can you please recommend me a proj object for projcrs function inside the projfwd function? I mean the code or wkt object (the projection), I don't know if a specific one would be more relevant than an other. My area is around 40 and 70°N, around UK. Thanks.
Sean de Wolski
Sean de Wolski 2021-8-10
编辑:Sean de Wolski 2021-8-10
The projection relies on where on the earth you are and the size of the area. For example, if working with Massachusetts data (MathWorks' headquarters location) I'd use the Massachusetts State Plane projection which is optimized for MA. Find it by web search literally: massachusetts state plane projection code - Bing
prj = projcrs(26986)
prj =
projcrs with properties: Name: "NAD83 / Massachusetts Mainland" GeographicCRS: [1×1 geocrs] ProjectionMethod: "Lambert Conic Conformal (2SP)" LengthUnit: "meter" ProjectionParameters: [1×1 map.crs.ProjectionParameters]
For Alaska, there are a bunch of them, so you'll have to pick one: Spatial Reference List -- Spatial Reference. Search as necessary for wherever in the world you are!
projcrs(3474)
ans =
projcrs with properties: Name: "NAD83(NSRS2007) / Alaska zone 7" GeographicCRS: [1×1 geocrs] ProjectionMethod: "Transverse Mercator" LengthUnit: "meter" ProjectionParameters: [1×1 map.crs.ProjectionParameters]

请先登录,再进行评论。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by