How to create a map mask using coordinates

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 ?

Answer 1

[MaHa]

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.

[Sean de Wolski]

this assumes that latitude and longitude are equal which they're not and so this will give wrong answers.

Answer 2

[Sean de Wolski]

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.)

[MaHa]

My end goal is to save data from my X matrix (same size as lat or lon) that is inside the polygon and play with it.

I am not sure I understand why the map projection is important, but thank you I'll give it a try next week !

[Sean de Wolski]

it's important because latitude and longitude are on the surface of a sphere, not on a flat cartesian coordinate system. Thus 1 latitude does not equal one longtitude and 1 latitude does not even equal 1 latitude if it's somewhere else. Thus you need to "flatten" it so that you can use cartesian things like poly2mask which assumes that pixels are square (a lat/lon cell is not square!)

[MaHa]

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]

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]