Finding the average value from polygon areas in a fishnet cell

Unless i missunderstand you this should work:

1. Intersect Fishnet and contours
2. Calculate a column of area-weighted values (value*shapearea)
3. Dissolve by fishnet ID and sum the area-weighted values
4. Calculate total averages by Dividing with shape areas

You can do this manually, with ModelBuilder or using the Python window and code below (change the four lines as indicated):

import arcpy
arcpy.env.overwriteOutput=1

fishnet = r'C:\Default.gdb\fishnet_4km' #Change to match your data
contours = r'C:\Default.gdb\contour' #Change to match your data
valuecolumn = 'AvgValue' #Change to match your data
outfc = r'C:\Default.gdb\fishnet_with_averages' #Change to match your data

#Intersect fishnet and contours
tempfc=r'in_memory\intersect'
arcpy.Intersect_analysis(in_features=[fishnet,contours], out_feature_class=tempfc)

#Calculate areaweighted values per intersection
arcpy.AddField_management(in_table=tempfc, field_name='AvgArea', field_type='DOUBLE')
with arcpy.da.UpdateCursor(tempfc,['AvgArea',valuecolumn,'SHAPE@AREA']) as cursor:
    for row in cursor:
        row[0]=row[1]*row[2]
        cursor.updateRow(row)

#Dissolve by fishnet ID and calculate sum of areasums
fishnetID = '{0}{1}'.format('FID_',arcpy.Describe(fishnet).name)
arcpy.Dissolve_management(in_features=tempfc, out_feature_class=outfc, 
                         dissolve_field=fishnetID, 
                         statistics_fields=[['AvgArea','SUM']])
arcpy.AddField_management(in_table=outfc, field_name='Areaweighted_average', field_type='DOUBLE')

#Calculate areaweighted average
with arcpy.da.UpdateCursor(outfc,['Areaweighted_average','SUM_AvgArea','SHAPE@AREA']) as cursor:
    for row in cursor:
        row[0]=row[1]/row[2]
        cursor.updateRow(row)

Inputs: Output:

The tool you need is Tabulate Intersection. You'll need to use a pivot table to transform the outputs to have one record per cell, and then join these results back to your fishnet feature class. The process is documented on the link above.

This can use a lot of memory, so make sure your files are local. A simple way to make it easier on your machine would be to split up the fishnet, and run it as a batch process. For more discussion and alternative python/R scripts see here

If I understand well, you should use the raster output of your spatial interpolation and zonal statistics if you want the most precise output. Don't be afraid with the discussion about exactness in Does IDW interpolation in ArcGIS Geostatistical analyst work as exact interpolation method? .

Statistically speaking, an exact interpolator is an interpolator where the interpolated value at the location of an input point is exactly the same as the value of this point. The answer then says that you are rarely exactly at the same location when you have pixels, therefore you will not have exactly the same value for a pixel and the "known" point that is inside this pixel (except in the center).

However, the continuous surface that is interpolated is always measured at the center of each pixel, which are then "exact" with respect to the interpolator. Pixels are therefore a good way to discretize an interpolated surface, and far more precise than polygons in your particular case (the hexagon is probably a circle that is simplified with 6 segments).