def zonal_statistics_for_grid(self, parameters, messages, dataset):

        zone_parameter = parameters[0]
        input_parameter = parameters[1]
        variable_parameter = parameters[2]
        output_parameter = parameters[3]
        type_parameter = parameters[4]
        ignore_parameter = parameters[5]
        zonefield_parameter = parameters[6]

        # Grab the parameter values
        dataset_name = input_parameter.valueAsText      # input netCDF filename
        output_filename = output_parameter.valueAsText  # output filename
        var_name = variable_parameter.valueAsText       # netCDF variable name

        # Extract variable from netCDF dataset
        arr_var = dataset.variable(var_name)

        # Check that spatial dimensions of input netCDF variable are regularly
        #   gridded and obtain cell size
        cell_size_x = None
        cell_size_y = None
        spat_list = list(dataset.space_dimension_names(var_name))
        assert len(spat_list) >= 2
        spat_var = dataset.variable(spat_list[0])
        blfirst = True
        if spat_var.size > 1:
            for i in range(spat_var.size - 1):
                cell_size_x = spat_var[i+1] - spat_var[i]
                if (not blfirst) and (cell_size_x != last_cell_size):
                    messages.addErrorMessage(
                        mds.messages.VARIABLE_NOT_REGULARLY_GRIDDED.format(
                            var_name, spat_list[0]))
                    raise arcpy.ExecuteError
                    break
                last_cell_size = cell_size_x
                blfirst = False
            cell_size_x = float(cell_size_x)
        spat_var = dataset.variable(spat_list[1])
        blfirst = True
        if spat_var.size > 1:
            for i in range(spat_var.size - 1):
                cell_size_y = spat_var[i+1] - spat_var[i]
                if (not blfirst) and (cell_size_y != last_cell_size):
                    messages.addErrorMessage(
                        mds.messages.VARIABLE_NOT_REGULARLY_GRIDDED.format(
                            var_name, spat_list[1]))
                    raise arcpy.ExecuteError
                    break
                last_cell_size = cell_size_y
                blfirst = False
            cell_size_y = float(cell_size_y)
        # If both x and y cell sizes are None, the dataset is a scalar.
        if (cell_size_x is None) and (cell_size_y is None):
            messages.addErrorMessage(
                mds.messages.VARIABLE_CONTAINS_INSUFFICIENT_VALUES.format(
                    var_name))
            raise arcpy.ExecuteError
        # If one of x and y is None, the dataset is a line, so we will assume
        #   the singleton dimension shares the other's cell size.
        if cell_size_x is None:
            cell_size_x = cell_size_y
            messages.addWarningMessage(
                ("Assuming variable dimension {} shares a cell size of {} " + \
                "with {}").format(spat_list[0], cell_size_y, spat_list[1]))
        if cell_size_y is None:
            cell_size_y = cell_size_x
            messages.addWarningMessage(
                ("Assuming variable dimension {} shares a cell size of {} " + \
                "with {}").format(spat_list[1], cell_size_x, spat_list[0]))

        # Figure out which spatial variable is x and which is y
        #   and then get corresponding number of columns and rows
        if dataset.convention.is_x_dimension_variable(dataset.variable(spat_list[0])):
            ncol = len(dataset.variable(spat_list[0]))
            nrow = len(dataset.variable(spat_list[1]))
            bltranspose = True
        else:
            ncol = len(dataset.variable(spat_list[1]))
            nrow = len(dataset.variable(spat_list[0]))
            cell_size_x, cell_size_y = cell_size_y, cell_size_x
            bltranspose = False

        # Grab extent from netCDF dataset
        #   i.e. [x_min, y_min, x_max, y_max]
        var_extent = dataset.extent(var_name)
        xl = float(var_extent[0])
        yl = float(var_extent[1])

        # Determine whether to flip the zone data over a given axis and
        #   transpose it to align with the netCDF data
        blflipx = False
        blflipy = False
        if cell_size_y > 0:
            blflipx = True
        if cell_size_x < 0:
            blflipy = True

        # Use absolute cell sizes
        cell_size_x = abs(cell_size_x)
        cell_size_y = abs(cell_size_y)

        # Convert extent to extent object for use later
        var_extent = arcpy.Extent(var_extent[0], var_extent[1], \
            var_extent[2], var_extent[3])

        # Was a temporary raster created?
        bltempRaster = False

        # Describe zone data
        desc = arcpy.Describe(zone_parameter.value)

        # Set environments
        arcpy.env.outputCoordinateSystem = desc.spatialReference
        temp_snap = arcpy.NumPyArrayToRaster(numpy.zeros((1,1)), \
            arcpy.Point(xl,yl), cell_size_x, cell_size_y)
        arcpy.env.snapRaster = temp_snap

        # Process input zone data
        valid_raster_types = ["RasterLayer", "RasterDataset"]
        valid_feature_types = ["ShapeFile","FeatureLayer", "FeatureDataset"]
        if desc.dataType in valid_feature_types:
            # Limitation of Feature to Raster
            if cell_size_x != cell_size_y:
                messages.addErrorMessage(
                    "The variable %s must use the same gridding " + \
                    "for both its spatial dimensions when feature zone " + \
                    "data is used." % var_name)
                raise arcpy.ExecuteError
            # Convert features
            temp_raster = 'in_memory\\temp'
            arcpy.FeatureToRaster_conversion(zone_parameter.value, \
                zonefield_parameter.value, temp_raster, cell_size_x)
            bltempRaster = True
            # NoData value
            noDataValue = arcpy.Describe(temp_raster).noDataValue
        elif desc.dataType in valid_raster_types:
            # Map non-Value field to Value if necessary and grab NoData value
            if zonefield_parameter.value != 'Value':
                in_raster = arcpy.sa.Lookup(zone_parameter.value, \
                    zonefield_parameter.value)
                # Raster created with Lookup will have default noDataValue
                # used for given .pixelType, but not assigned to .noDataValue
                if '8' in in_raster.pixelType:
                    noDataValue = -128  # S8
                else:
                    noDataValue = -2147483647 # S32
            else:
                in_raster = zone_parameter.valueAsText
                noDataValue = desc.noDataValue
            # Resample if necessary
            if (desc.meanCellHeight != cell_size_y) or \
                (desc.meanCellWidth != cell_size_x):
                temp_raster = 'in_memory\\temp'
                arcpy.Resample_management(in_raster, temp_raster,
                    "%f %f" % (cell_size_x, cell_size_y))
                bltempRaster = True
            else:
                temp_raster = in_raster
        else:
            messages.addErrorMessage("Not a valid zone dataset.")
            raise arcpy.ExecuteError

        # Make raster object from zone data if not already
        if not isinstance(temp_raster, arcpy.Raster):
            temp_raster = arcpy.Raster(temp_raster)

        # Limit to lon [-180, 180] and lat [-90, 90]
        def convert_coord(num, islon = True):
            lim = 180 if islon else 90
            if num > lim:
                out = num % -lim
            elif num < -lim:
                out = num % lim
            else:
                out = num
            if not islon:
                mult = -1.0 if num < 0 else 1.0
                out = mult * abs(out)
            return out

        # Determine whether input zone data has a PCS or GCS
        if (temp_raster.spatialReference.exporttostring() == \
                temp_raster.spatialReference.GCS.exporttostring()):
            # Test for overlap between zone and netCDF data
            if temp_raster.extent.disjoint(var_extent):
                # Try with coordinates in lon[-180,180] and lat[-90,90] ranges
                var_extent_temp = arcpy.Extent(
                    convert_coord(var_extent.XMin, True),
                    convert_coord(var_extent.YMin, True),
                    convert_coord(var_extent.XMax, False),
                    convert_coord(var_extent.YMax, False))
                xl = var_extent_temp.XMin
                yl = var_extent_temp.YMin
                # Check again and throw warning on failure
                if temp_raster.extent.disjoint(var_extent_temp):
                    messages.addWarningMessage("Data extents don't overlap.")
        else:
            if (not temp_raster.extent.overlaps(var_extent)) and (
                    not temp_raster.extent.equals(var_extent)):
                messages.addWarningMessage("Data extents don't overlap.")

        # Extract zone data from raster
        arr_zone = arcpy.RasterToNumPyArray(temp_raster, arcpy.Point(xl,yl), \
            ncol, nrow)
        # Flip zone raster to accord with how data is stored
        if blflipx:
            arr_zone = numpy.flipud(arr_zone)
        if blflipy:
            arr_zone = numpy.fliplr(arr_zone)
        if bltranspose:
            arr_zone = numpy.transpose(arr_zone)
        # Delete temporary raster if it exists
        if bltempRaster:
            arcpy.Delete_management(temp_raster)

        # Write new dataset
        with netCDF4.Dataset(output_filename, mode="w", \
                format=dataset.dataset.file_format) as newdataset:

            # Copy global attributes
            for attribute_name in dataset.dataset.ncattrs():
                newdataset.setncattr(attribute_name, dataset.dataset.getncattr(
                    attribute_name))
            newdataset.Source_Software = "Esri ArcGIS"
            # Loop over dimensions in original dataset
            for dim_name in dataset.dimension_names():
                # Create dimension
                old_dim = dataset.dimension(dim_name)
                newdataset.createDimension(dim_name, \
                    len(old_dim) if not old_dim.isunlimited() else None)
                if dim_name in dataset.variable_names():
                    # Create variable corresponding to dimension
                    dim_var = dataset.variable(dim_name)
                    new_var = newdataset.createVariable(dim_name, dim_var.dtype, \
                        dim_name)
                    # Copy attributes
                    for attribute_name in dim_var.ncattrs():
                        new_var.setncattr(attribute_name, dim_var.getncattr(
                            attribute_name))
                    # Copy data
                    new_var[:] = dim_var[:]
            # Loop over any dependent variables
            dim_names = []
            if 'coordinates' in arr_var.ncattrs():
                dim_names = str(arr_var.getncattr('coordinates')).split()
            for dim_name in dim_names:
                if dim_name not in dataset.dimension_names():
                    # Create variable corresponding to dimension
                    dim_var = dataset.variable(dim_name)
                    new_var = newdataset.createVariable(dim_name, dim_var.dtype, \
                        dim_var.dimensions)
                    # Copy attributes
                    for attribute_name in dim_var.ncattrs():
                        new_var.setncattr(attribute_name, dim_var.getncattr(
                            attribute_name))
                    # Copy data
                    new_var[:] = dim_var[:]
            # Check if variable is packed
            ispacked = hasattr(arr_var,'scale_factor') or hasattr(arr_var, 'add_offset')
            # Determine whether given statistic is categorical or not
            new_var_dtype = False
            for stat in self.statistics:
                if type_parameter.valueAsText in stat[1].keys():
                    new_var_dtype = stat[2]
                    break
            type_changed = False
            if (('int' in str(arr_var.dtype)) and (not new_var_dtype)) and (
                    not ispacked):
                type_changed = True
            # Create new variable
            #   retaining type of original variable unless an integer variable is
            #   selected with a floating-point (non-categorical) statistic
            new_var = newdataset.createVariable(var_name, 'f8' if type_changed \
                else arr_var.dtype, dataset.variable_dimension_names(var_name))
            # Copy attributes
            attr_names = ['_FillValue', 'missing_value', 'valid_min', 'valid_max', \
                'valid_range']
            for attribute_name in arr_var.ncattrs():
                # Handle _FillValue type mismatch error when changing dtype
                # of variable; also change for related (but non-error-causing)
                # attributes
                if type_changed and (str(attribute_name) in attr_names):
                    new_var.setncattr(attribute_name, \
                        arr_var.getncattr(attribute_name).astype('f8'))
                else:
                    new_var.setncattr(attribute_name, arr_var.getncattr(
                        attribute_name))

            # Slice generator
            # Adapted from: http://code.activestate.com/recipes/
            #   502194-a-generator-for-an-arbitrary-number-of-for-loops/
            # Originally written by Steven Bethard
            def multi_for(iterables):
                if not iterables:
                    yield ()
                else:
                    if isinstance(iterables[0], collections.Iterable):
                        for item in iterables[0]:
                            for rest_tuple in multi_for(iterables[1:]):
                                yield (item,) + rest_tuple
                    else:
                        for rest_tuple in multi_for(iterables[1:]):
                            yield (iterables[0],) + rest_tuple

            # Generate slice list
            dim_shape = arr_var.shape
            dim_name = arr_var.dimensions
            dim_slices = []
            dim_prod = numpy.int(1)
            for dim_index, dim_item in enumerate(dim_name):
                if not dim_item in spat_list:
                    dim_prod *= dim_shape[dim_index]
                    dim_slices.append(xrange(dim_shape[dim_index]))
                else:
                    dim_slices.append(Ellipsis)

            # Get zone values
            if noDataValue is not None:
                zones = numpy.setdiff1d(numpy.unique(arr_zone), [noDataValue])
            else:
                zones = numpy.unique(arr_zone)

            # Progress bar
            dim_prod *= zones.size
            arcpy.ResetProgressor()
            arcpy.SetProgressor('step', 'Calculating...', 0, dim_prod, 1)

            # Loop through slices
            for varslice in multi_for(dim_slices):

                # Take slice from data
                if type_changed:
                    arr_slice = arr_var[varslice].astype('f8')
                else:
                    arr_slice = arr_var[varslice]

                # Allocate slice to save result
                arr_newslice = numpy.ma.masked_all_like(arr_slice)

                # Create masked array
                if hasattr(arr_slice, 'mask'):
                    arr_mask = arr_slice.mask.copy()
                else:
                    arr_mask = numpy.ma.make_mask_none(arr_slice.shape)
                    arr_slice = numpy.ma.asarray(arr_slice)

                # Loop over zones
                for z in zones:
                    # Perform for zone when NoData is ignored
                    #   or when there is no NoData within zone
                    if ignore_parameter.value or (not numpy.any((arr_zone==z) & arr_mask)):
                        # Reset mask
                        arr_slice.mask = arr_mask
                        # Apply zone mask
                        arr_slice[arr_zone!=z] = numpy.ma.masked
                        # Apply operator to zone masked slice
                        result = self.calculate_statistic(arr_slice[:], \
                            type_parameter.valueAsText)
                        # Fill zone with result
                        arr_newslice[(arr_zone==z)] = result

                    arcpy.SetProgressorPosition()

                # Write slice
                new_var[varslice] = arr_newslice

        arcpy.ResetProgressor()

        return
Compare with another reference
dataset/ETHPy150Open Esri/solutions-geoprocessing-toolbox/suitability/toolboxes/scripts/MultidimensionSupplementalTools/MultidimensionSupplementalTools/Scripts/mds/tools/multidimensional_zonal_statistics.py/MultidimensionalZonalStatistics.zonal_statistics_for_grid--crossover with--dataset/ETHPy150Open AppScale/appscale/AppServer/lib/django-1.5/django/dispatch/saferef.py/BoundMethodWeakref.__eq__