map_coverage_gs.asv

% function for showing map coverage and ground stations locaiton
% some bugs to be fixed

function [] = map_coverage_gs()
    
    latlim = [22.00417 25.12825];
    lonlim = [118.31833 121.753];
    spacingInLatLon = 1; % 基站間距離，越小覆蓋率越高但計算時間較長(單位degrees)
    EPSG_code = 3826; % (TWD97-EPSG:3826/3830)
    proj = projcrs(EPSG_code); % projected coordinate reference system
    spacingInXY = deg2km(spacingInLatLon)*1000; % 相對於spacingInLatLon轉換為xy座標系統(單位:meters)
    [xlim,ylim] = projfwd(proj,latlim,lonlim);
    R = maprefpostings(xlim,ylim,spacingInXY,spacingInXY);
    [X,Y] = worldGrid(R);
    [gridlat,gridlon] = projinv(proj,X,Y);
    
    TW_landareas = readgeotable("E:\MATLAB\碩士論文\Coverage_Maps_for_Satellite_Constellation\gadm36_TWN_shapread\gadm36_TWN_0.shp");
    % ROC = landareas(landareas.Name == "Taiwan",:);
    
    T = geotable2table(TW_landareas,["Latitude","Longitude"]);
    [landlat,landlon] = polyjoin(T.Latitude,T.Longitude);
    
    bufwidth = 1;
    [landlatb,landlonb] = bufferm(landlat,landlon,bufwidth,"outPlusInterior");
    TW_liab = geopolyshape(landlatb,landlonb);
    
    gridpts = geopointshape(gridlat,gridlon);
    inregion = isinterior(TW_liab,gridpts);
    gslat = gridlat(inregion);
    gslon = gridlon(inregion);
    
    gs = groundStation(sc,gslat,gslon);
    
    % Add Ground Station Receivers (地面接收站獲取訊號後透過基地台轉傳)
    isotropic = arrayConfig(Size=[1 1]);
    tx = receiver(gs,Antenna=isotropic);
    
    %% have some problem to be fix
    pattern(tx,Size=500000);
    
   

end