% This code finds all the SAR swaths within a given ROI
clear; close all;

% input ROI : either min & max lat, lon, or an input shapefile
% minlon=145; minlat=-20; maxlon=175; maxlat=17;
% polyROI = polyshape([minlon maxlon maxlon minlon],[minlat minlat maxlat maxlat]); 
shpROIname = 'F:\Projects\2023_NP_mesoscale_toposynth\shapefiles\ROI_NP_mesoscale_proj.shp';
shpROI = shaperead(shpROIname);
polyROI = polyshape(shpROI.X,shpROI.Y);

% plot the ROI just for fun. 
figure; hold on; plot(polyROI)

% SAR swath outlines (this includes all SAR swaths, from PDS)
pathdata = 'F:\Titan_data\RADAR\footprints\2021_PDS\';
shpswaths = shaperead([pathdata 'swath_outlines.shp']);
nswaths = length(shpswaths);

% now find the swaths intersecting the ROI
swaths_in=shpswaths(1); kk=0;
for ss = 1:nswaths
    outline_ii = shpswaths(ss);
    polyshp = polyshape(outline_ii.X,outline_ii.Y); 
    polyout = intersect(polyROI,polyshp);
    if polyout(1).NumRegions>0
        kk=kk+1; 
        swaths_in(kk)=shpswaths(ss);
        % plot them for fun
        plot(polyshp); 
    end
end
