Alexandre, boa tarde!
Sugiro que você trabalhe com o mapa de bacias "dissolvidas" para o nível 2. Para isso você pode fazer uso de maptools::unionSpatialPolygons(). Para fins de visualização você pode ainda utilizar outras funções de simplificação.
Feito isso, além de visualizar mais rápido, você poderá obter o centróide dos polígonos dissolvidos pra locar os gráficos de torta. Depois de rodar o sp::over() você precisa tabular as frequências por bacia e adicionar a coordenada de cada uma.
A entrada do add.pie() é individual, então você poderá pensar num laço ou uma função da família apply pra entradas múltiplas.
Atte.,