# How to create a map mask using coordinates

39 views (last 30 days)
MaHa on 6 Aug 2021
Edited: Sean de Wolski on 10 Aug 2021
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:
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 on 9 Aug 2021
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
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.
Sean de Wolski on 9 Aug 2021
this assumes that latitude and longitude are equal which they're not and so this will give wrong answers.

Sean de Wolski on 6 Aug 2021
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.)
Sean de Wolski on 10 Aug 2021
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]