[GIS] Spatial join of multiple different point layers keeping all attributes in ArcGIS

Question

I have 28 fairly large point data sets (1/6 degree global grid), all of different extents but with the same point location (previously extracted as centroids from raster pixel).

Each point layer has a unique value (grid_code in example below) per point but as the spatial extents for each point file vary, FID and pointid do not match, i.e. there is no common value to join on.

Using ArcGIS Desktop (10.5) I'd like to combine all 28 files to one by location as the Spatial Join tool does (a reverse one-to-many), but create one file that contains the full spatial extent of all files combined and the respective values attributes (preferably with the origin layer names) per point:

I have tried all varieties of merge, join, relate and spatial join but they either do not combine the attributes, do not join by location or there is no possibility to do it as a batch function.

Answer 1

You can do this using some Python code and the pandas module which is included in 10.5. All input Point feature classes need to be in one database (and no other features). They all need to have a field called grid_code. A field will be added to each of the feature classes named gcode+feature class name, for example gcode_Point1. This is then given the value of grid_code. Then all fcs are merged together, exported to pandas, grouped by xy coordinates (which need to be exactly the same for overlapping points for it to work. But this can be adjusted by rounding the coordinates) so overlapping Points are combined and all non-overlapping Points kept. This is then exported to a new feature class.

Backup your data Before you try it! The code can be executed in the Python window after you change Point_database.

import arcpy
import pandas as pd

point_database = r'C:\Test.gdb' #All input feature classes here, nothing else. Change to match the name of your database
output_fc = 'MegaMerge'
arcpy.env.workspace = point_database

all_features = arcpy.ListFeatureClasses()

for feature in all_features:
    fname = 'gcode_{0}'.format(feature)
    arcpy.AddField_management(in_table=feature, field_name=fname, field_type='DOUBLE')
    with arcpy.da.UpdateCursor(feature,['grid_code',fname]) as cursor:
        for row in cursor:
            row[1]=row[0]
            cursor.updateRow(row)

arcpy.Merge_management(inputs=all_features, output='all_points')
fieldlist = [f.name for f in arcpy.ListFields('all_points') if f.name.startswith('gcode')]
df = pd.DataFrame.from_records(data=arcpy.da.SearchCursor('all_points',['SHAPE@XY']+fieldlist), columns=['SHAPE@XY']+fieldlist)
df_grouped = df.groupby('SHAPE@XY').sum()

arcpy.CreateFeatureclass_management(out_path=arcpy.env.workspace, out_name=output_fc, 
                                   geometry_type='POINT', 
                                   spatial_reference=arcpy.Describe(all_features[0]).spatialReference)

for name in list(df_grouped.columns.values):
    arcpy.AddField_management(in_table=output_fc, field_name=name, field_type='DOUBLE')

icur = arcpy.da.InsertCursor(output_fc, ['SHAPE@XY']+fieldlist)
for row in df_grouped.itertuples(index=True, name=None):
    icur.insertRow(row)
del icur