Discussion :

[thecrazydonut]

Bonjour à tous,

Je cherche à identifier les points d'une grille qui chevauchent un polygone défini dans un shapefile. Plus précisément, je voudrais pouvoir retourner les indices des vecteurs contenant les longitudes et les latitudes qui sont compris à l'intérieur de ce polygone. Ou encore mieux, pouvoir définir une distance à partir de laquelle un point est considéré comme à l'intérieur du polygone.

J'ai produit le shapefile à l'aide d'arcGIS. Il contient un polygone "déprojeté", c'est à dire que la position des points est exprimée en degrés. Les grilles sont issues de centres de prévisions météo et stockées dans des fichiers NetCDF. Il peut s'agir de grilles régulières (cas de ECMWF dans le code dessous), ou correspondant à une projection (projection Lambert pour NOAA).

Une petite figure pour illustrer le problème.
La ligne bleue représente les limites de mon polygone, les points noirs les nœuds de la grille ECMWF, les points rouge les noeuds de la grille NOAA.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
% Import shape file
[S,A]=shaperead('./023402_shp/023402.shp');
 
% Import NetCDF NOAA
ncid = netcdf.open('./NetCDF_files/noaa.nc','NOWRITE');
lat0noaa = netcdf.getVar(ncid, netcdf.inqVarID(ncid,'g3_lat_0'));
lon0noaa = netcdf.getVar(ncid, netcdf.inqVarID(ncid,'g3_lon_1'));
netcdf.close(ncid)
 
% Import NetCDF ECMWF
ncid = netcdf.open('./NetCDF_files/ecmwf.nc','NOWRITE');
lat0ecmwf = netcdf.getVar(ncid, netcdf.inqVarID(ncid,'lat_0'));
lon0ecmwf = netcdf.getVar(ncid, netcdf.inqVarID(ncid,'lon_0'));
netcdf.close(ncid)
 
% Plots
plot(S.X,S.Y,'b'); hold on
plot(lon0noaa,lat0noaa,'LineStyle','none','Marker','.','Color','r')
plot(repmat(lon0ecmwf-360, [1, size(lat0ecmwf,1)])',...
    repmat(lat0ecmwf, [1, size(lon0ecmwf,1)]),...
    'LineStyle','none','Marker','.','Color','k')
set(gca,'XLim',[-72 -70],'YLim',[45 47])

Je mets en pièce jointe des données pour faire tourner le code (seulement la grille rouge et le shapefile, sinon la pièce jointe est trop lourde).

Merci pour vos suggestions !

[Jerome Briot]

Regarde la fonction inpolygon ou ses alternatives sur le File Exchange

[le fab]

salut

as tu essayer la fonction inpolygon ?

Fabien

[Jerome Briot]

J'ai oublié de préciser que dans ce genre de problème, tu peux initialement réduire le domaine de recherche en excluant les points qui sont en dehors du rectangle englobant le polygone :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
idx = x<min(S.X(:)) | x>max(S.X(:)) | y<min(S.Y(:)) | y>max(S.Y(:));
 
x(idx) = [];
y(idx) = [];

Tu peux aussi remplacer le rectangle par le polygone convexe englobant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
k = convhull(xpolygon, ypolygon);
 
in = inpolygon(x, y, x(k), y(k);
 
x(~in) = [];
y(~in) = [];

C'est l'idée…

[thecrazydonut]

Merci, je suis passé à côté de cette fonction lors de mes recherches ! C'était ça qu'il me fallait.
Par contre pas d'option pour spécifier une distance à partir de laquelle un point est "attrapé"

[le fab]

bein en fait soit le point est dedans, soit il ne l'est pas
du coup pour attraper des points qui ne son pas dedans mis qui sont pas loin, il faut tricher
par exemple en faisant l'union de plusieurs inpolygone ou l'on décale les grilles X,Y dans les 4 directions