matlab - 根据经纬度获取大陆

Question

<分区>

Answer 1

我完成这段代码只是为了获得大陆。它不是 100% 准确，但对大多数地区来说都很好。

我这样做是为了按照 PVsyst 的要求识别大陆，因此我区分了澳大利亚和太平洋地区，但这很容易改变:

如果您不想区分澳大利亚和太平洋地区，只需删除澳大利亚部分(LatAus、LonAus，将其从 contientIndex 和大陆列表中删除)。

function Continent = GetContinent(Lat,Lon,plotTrueFalse)

LatNAm = [90       90  78.13   57.5  15  15  1.25  1.25  51  60    60];
LonNAm = [-168.75 -10 -10     -37.5 -30 -75 -82.5 -105  -180 -180 -168.75];
LatNA2 = [51    51  60];
LonNA2 = [166.6 180 180];

LatSAm = [1.25   1.25  15  15 -60 -60];
LonSAm = [-105 -82.5  -75 -30 -30 -105];

LatEur = [90   90  42.5 42.5 40.79 41 40.55 40.40 40.05 39.17 35.46 33   38  35.42 28.25 15  57.5  78.13];
LonEur = [-10 77.5 48.8 30   28.81 29 27.31 26.75 26.36 25.19 27.91 27.5 10 -10  -13   -30 -37.5 -10];

LatAfr = [15  28.25 35.42 38 33   31.74 29.54 27.78 11.3 12.5 -60 -60];
LonAfr = [-30 -13   -10 10 27.5 34.58 34.92 34.46 44.3 52    75 -30];

LatAus = [-11.88 -10.27 -10 -30    -52.5 -31.88];
LonAus = [110      140  145 161.25 142.5  110];

LatAsi = [90   42.5 42.5 40.79 41 40.55 40.4  40.05 39.17 35.46 33   31.74 29.54 27.78 11.3 12.5 -60 -60 -31.88 -11.88 -10.27 33.13 51    60  90];
LonAsi = [77.5 48.8 30   28.81 29 27.31 26.75 26.36 25.19 27.91 27.5 34.58 34.92 34.46 44.3 52   75  110  110   110    140    140   166.6 180 180];
LatAs2 = [90    90      60      60];
LonAs2 = [-180 -168.75 -168.75 -180];

LatAnt = [-60 -60 -90 -90];
LonAnt = [-180 180 180 -180];

inNAm         = inpolygon(Lat,Lon,LatNAm,LonNAm) + inpolygon(Lat,Lon,LatNA2,LonNA2);
inEur         = inpolygon(Lat,Lon,LatEur,LonEur);
[inSAm,onSAm] = inpolygon(Lat,Lon,LatSAm,LonSAm);
[inAfr,onAfr] = inpolygon(Lat,Lon,LatAfr,LonAfr);
[inAus,onAus] = inpolygon(Lat,Lon,LatAus,LonAus);
[inAsi,onAsi] = inpolygon(Lat,Lon,LatAsi,LonAsi);
[inAs2,onAs2] = inpolygon(Lat,Lon,LatAs2,LonAs2);

ContinentIndex  = 1 +(inAfr||onAfr)*1 + (inAus||onAus)*2 + (inAsi||onAsi||inAs2||onAs2)*3 + inEur*4 + inNAm*5 + (inSAm||onSAm)*6;
ContinentList   = {'Pacific','Africa','Australia','Asia','Europe','North_America','South_America'};

Continent       = ContinentList{ContinentIndex};
%% If you want to show the figure, set plotTrueFalse to true or to ~=0.
if plotTrueFalse
    geoshow('landareas.shp', 'FaceColor', [1 1 1]);
    axis([-180 180 -90 90])
    set(gca,'xtick',-180:15:180,'ytick',-90:15:90)
    grid on;
    hold on;
    patch(LonNAm,LatNAm,'r');
    patch(LonNA2,LatNA2,'r');
    patch(LonSAm,LatSAm,'g');
    patch(LonEur,LatEur,'b');
    patch(LonAfr,LatAfr,'k');
    patch(LonAus,LatAus,[1 0.5 0]);
    patch(LonAsi,LatAsi,'y');
    patch(LonAs2,LatAs2,'y');
    patch(LonAnt,LatAnt,[0 1 1]);
    alpha(0.5);
end