function shps = loadShapeFile(filename,boundary,params)

is_geodetic = params.is_shp_geodetic;

% Read data
shps = shaperead(filename);

%Split by NAN and convert to UTM
nshps = numel(shps);
new_shps=[];
for s=1:nshps
    orig_shp = shps(s);
    if(is_geodetic)
        [orig_shp.X,orig_shp.Y,orig_shp.Zone] = ll2utm(orig_shp.Y,orig_shp.X);
    end
    
    %Ignore vectors that are completely outside boundary
    if(~isempty(boundary))
        inside = isPointInPolygon([orig_shp.X',orig_shp.Y'],boundary);
        if(~any(inside))
            continue;
        end
    end
    
    if(~isnan(orig_shp.X(end)))
        orig_shp.X(end+1) = nan;
        orig_shp.Y(end+1) = nan;
    end
    
    splits = find(isnan(orig_shp.X));
    nsplits = numel(splits);
    splits = [0,splits];
    for k=1:nsplits
        shp = orig_shp;
        
        first = splits(k)+1;
        last = splits(k+1)-1;
        
        shp.X = shp.X(first:last);
        shp.Y = shp.Y(first:last);
        
        shp.BoundingBox(1) = min(shp.X);
        shp.BoundingBox(2) = min(shp.Y);
        shp.BoundingBox(3) = max(shp.X);
        shp.BoundingBox(4) = max(shp.Y);
        
        shp.offsets = zeros(1,2);
        
        new_shps = [new_shps;shp];
    end
end
shps = new_shps;

end


