How to create a map mask using coordinates

21 views (last 30 days)
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 ?

Accepted Answer

MaHa
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
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 Comment
Sean de Wolski
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.

Sign in to comment.

More Answers (1)

Sean de Wolski
Sean de Wolski on 6 Aug 2021
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 Comments
MaHa
MaHa on 10 Aug 2021
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 on 10 Aug 2021
Edited: 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]

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!