I have three polygon shapefiles that each represent water (lakes, ponds, rivers) for three different time periods (1990, 2000, 2010). Each shapefile has over 1 million polygons/attributes and were generated from rasters, using a pretty long process in Arcpy (rough summary of the process: region grouping using all 8 neighbors, raster to polygon conversion, dissolve along the gridcode created by region grouping). 

What I want to do is see how the water bodies have changed across these time periods, in particular, seeing how their areas have changed. In some cases, a lake/pond that existed in 1990 (for example) could have partially dried up and turned into 2 lake/ponds, or likewise, two lakes/ponds from 1990 could now be one polygon. The output doesn't have to be a shapefile, I just would (eventually) like a list of all the polygons represented in all of the shapefiles, and their corresponding areas for each time period.

So far, what I tried to do is a Spatial Join between the 1990 and 2000 shapefiles, then using that result (90_00_join.shp) in another Spatial Join with the 2010 shapefile (result is 90_00_10_join.shp). I wanted to capture the cases I mentioned above, so I used the Join one to many option. I also want all polygons from each time period to be present in the final output. So if a polygon (a.k.a. lakes/ponds) that existed in 1990 and/or 2000 but doesn't exist in 2010, I still want it in the final output. If a polygon didn't exist in earlier years but does in 2010, I also want it in the final output.

However, Spatial Join only allows for me to "Keep All Target Features." As a result, my final Spatial Join output has the same number of attributes as the 1990 input shapefile (the one that was used as the Target Feature in making 90_00_join.shp), although there are about 60,000 more attributes in the 2000 shapefile. So clearly these attributes from 2000 are getting left behind, I'm assuming because they are the 'Join feature' and not the Target feature. The same then happens after I join 90_00_join.shp with the 2010 shapefile.


Is there a way to get what I want using Spatial Join? If not, what about in python/arcpy by looping through the features and getting the area of each polygon and the polygons from the other shapefiles that intersect each feature? If not in Arc, what about QGIS or some other open source GIS software? The thing to consider is the number of polygons each shapefile has (about 1.2 - 1.3 million polygons).