Comparing areas of polygon over time

Question

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).

EDIT: I should have mentioned that the problem I'm having is that the overlapping polygons do not have the same ID. So there is no way to compare the polygons across the years by just looking at the attribute tables. What I would like is to eventually have a table with ALL the polygons that existed throughout each time period, and their areas for each of the time period. This could normally be achieved using a spatial join, however what I can't figure out is how to get the output to include all polygons, not just those from the Target feature.

Answer 1

The biggest problem seems like the fact that you don't have a primary key to relate the three datasets to one another. Try appending all the polygons from each dataset into a single feature class and dissolve the result into single part polygons. This result would represent the combined boundary of each lake over time. On this layer (ID layer), add a field (like "Lake_ID", type LONG) and field calculate an incremental ID (field calculator, answer to this question: Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator?).

After, for each of your inputs, add your Lake_ID field, spatial join on the ID layer, then use field calculator to insert the Lake_ID value from the ID layer into the same field in your input. Once this is done, Lake_ID for your input layers will be a secondary key. You can then dissolve on this field, and use the "Statistics Fields" option to calculate the sums of each of your value. For example, if you've already calculated the area (or partCount using !SHAPE!.partCount) of each individual polygon, you could SUM this field to get the total in your output.

At this point, each dataset representing the three different years should have the same number of rows. You can then perform a table join (one to one) on all of them, and export the join to a new table. Alternatively, you could use Make Query Table (http://pro.arcgis.com/en/pro-app/tool-reference/data-management/make-query-table.htm) to join the three datasets, and choose which fields you want in your output. If you do this, you'll have to export the dataset to a shapefile/feature class to keep the results. With the Make Query table tool, if you keep a SHAPE attribute from one of the inputs, the result will be a query layer with geometries.

Answer 2

Re class 1st raster to 1,0. Second 2,0. Third 4,0.

Use raster calculator to find total.

Use 1st year to tabulate area of totals. For each polygon you'll get areas of for values in range 0, 1, 3, 5, 7. I hope meaning is obvious.

Same can be done for other years.

This is useful output unless you are simply talking about merging them