def zonal_statistics_for_discrete(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

        # Get associated station, x, and y variables
        (stat_variable, x_variable, y_variable) = self.get_dependent_variables(
            dataset, var_name)
        if (stat_variable is None) or (x_variable is None) or (y_variable is None):
            messages.addErrorMessage('%s is not a station variable.' % var_name)
            raise arcpy.ExecuteError

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

        # Extract coordinates attribute
        dim_names = []
        if 'coordinates' in arr_var.ncattrs():
            dim_names = str(arr_var.getncattr('coordinates')).split()

        # Get coordinate arrays
        X = dataset.variable(x_variable)
        Y = dataset.variable(y_variable)

        # Make array of station indices
        if stat_variable in dataset.variable_names():
            arr_stat = dataset.variable(stat_variable)[:]
        else:
            arr_stat = numpy.arange(len(dataset.dimension(stat_variable)))

        # Empty dictionary to hold stations falling within each zone
        zone_stations = {}
        #   where the keys are values of the zone field and the key values
        #   are each a list of station indices contained within that zone

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

        # Was a temporary feature class created?
        bltemppoly = False

        # Convert raster to polygon if necessary
        valid_raster_types = ["RasterLayer", "RasterDataset"]
        valid_feature_types = ["ShapeFile","FeatureLayer", "FeatureDataset"]
        if desc.dataType in valid_raster_types:
            temp_poly = "in_memory\\temp"
            arcpy.RasterToPolygon_conversion(in_raster, temp_poly, \
                "NO_SIMPLIFY", zonefield_parameter.value)
            temp_field = "GRIDCODE"
            bltemppoly = True
        elif desc.dataType in valid_feature_types:
            temp_poly = zone_parameter.value
            temp_field = zonefield_parameter.value
        else:
            messages.addErrorMessage("Not a valid zone dataset.")
            raise arcpy.ExecuteError

        num_recs = int(arcpy.GetCount_management(temp_poly).getOutput(0))
        num_recs *= numpy.prod(arr_stat.shape)

        arcpy.SetProgressor("step", "Building zone feature look-up...", 0, num_recs, 1)

        # Find zones associated with each feature
        desc = arcpy.Describe(temp_poly)
        rows = arcpy.SearchCursor(temp_poly)
        for row in rows:
            feat = row.getValue(desc.ShapeFieldName)
            zone = row.getValue(temp_field)
            for station in arr_stat:
                if feat.contains(arcpy.Point(float(X[station]),float(Y[station]))):
                    if not zone in zone_stations.keys():
                        zone_stations[zone] = []
                    zone_stations[zone].append(station)
                arcpy.SetProgressorPosition()
        del rows, row
        arcpy.ResetProgressor()

        # Delete temporary feature class if it exists
        if bltemppoly:
            arcpy.Delete_management(temp_poly)

        # 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 dependent variables in original dataset
            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 dim_item != stat_variable:
                    dim_prod *= dim_shape[dim_index]
                    dim_slices.append(xrange(dim_shape[dim_index]))

            # Get position of statistic
            stat_position = dim_name.index(stat_variable)

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

            # Loop through zones
            for zone in zone_stations:

                # Grab station indices associated with zone
                zone_indices = zone_stations[zone]

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

                    # Insert station indices into slice
                    varslice = list(varslice)
                    varslice.insert(stat_position, zone_indices)

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

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

                    # Perform for zone when NoData is ignored
                    #   or when there is no NoData within zone
                    if ignore_parameter.value or (not numpy.any(arr_mask)):

                        # Apply mask
                        arr_slice.mask = arr_mask

                        # Apply operator to masked slice
                        result = self.calculate_statistic(arr_slice[:], \
                            type_parameter.valueAsText)

                        # Write slice if possible
                        if not isinstance(result, numpy.ma.core.MaskedConstant):
                            new_var[varslice] = result

                    arcpy.SetProgressorPosition()

        arcpy.ResetProgressor()

        return
auto generate default config file, and stage it into the profile.
dataset/ETHPy150Open Esri/solutions-geoprocessing-toolbox/suitability/toolboxes/scripts/MultidimensionSupplementalTools/MultidimensionSupplementalTools/Scripts/mds/tools/multidimensional_zonal_statistics.py/MultidimensionalZonalStatistics.zonal_statistics_for_discrete--crossover with--dataset/ETHPy150Open ipython/ipython-py3k/IPython/core/application.py/BaseIPythonApplication.stage_default_config_file