def zonal_statistics_as_table_for_discrete(self, parameters, messages, dataset):
        """Performs zonal statistics as table assuming dataset is a CF 1.6
        compliant discrete sampling geometry dataset of 'point', 'timeseries',
        or 'profile' feature type using an orthogonal or incomplete
        multidimensional array data representation."""

        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)

        # 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 stat[1].keys() in type_parameter.valueAsText:
                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

        # 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_dims = []
        dim_prod = 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]))
                dim_dims.append((str(dim_item), numpy.int))

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

        blstructcreated = False

        def hstack_struct(arr1, arr2):
            # Will only work correctly if arr1 and arr2 only have one record
            dt = numpy.lib.recfunctions.zip_descr((arr1, arr2), flatten=True)
            return numpy.array(arr1.tolist() + arr2.tolist(), dtype = dt)

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

        # Loop through zones
        index = 0
        valid_indices = []
        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 zone masked slice
                    result = self.calculate_statistic(arr_slice[:], \
                        type_parameter.valueAsText)

                    # Structure results into row
                    struct_data = []
                    for it_ind, it in enumerate(list(varslice)):
                        if it_ind != stat_position:
                            if dim_name[it_ind] in dataset.variable_names():
                                # Dimension with variable
                                struct_data.append(dataset.variable(dim_name[it_ind])[it])
                            else:
                                # Dimension without variable
                                struct_data.append(it)
                    struct_data.append(zone)
                    struct_data.append(numpy.ma.count(arr_slice))
                    if numpy.isscalar(result):
                        # Single-valued result
                        struct_data.append(result)
                    else:
                        # Multivalued result
                        struct_data += [it2 for it in result.tolist() for it2 in it]

                    # Create table when necessary
                    if not blstructcreated:
                        struct_dtype = dim_dims[:]
                        struct_dtype.append((str(zonefield_parameter.value),numpy.int))
                        struct_dtype.append(('COUNT', numpy.int))
                        if numpy.isscalar(result):
                            # Single-valued result
                            struct_dtype.append((str(type_parameter.valueAsText), \
                                numpy.float))
                        else:
                            # Multivalued result
                            struct_dtype += [(it, numpy.float) for it in list(
                                result.dtype.names)]
                        struct = numpy.zeros((dim_prod,), numpy.dtype(struct_dtype))
                        blstructcreated = True

                    # Insert row
                    struct[index] = numpy.array(tuple(struct_data), \
                        numpy.dtype(struct_dtype))
                    valid_indices.append(index)
                #else:
                #    # Passing drops zones with NoData in them from table
                #    pass
                index += 1
                arcpy.SetProgressorPosition()

        arcpy.ResetProgressor()

        if 'time' in dim_name:
            # Output zonal statistics to temporary table
            temp_table = os.path.join(tempfile.gettempdir(), 'temp_table.dbf')
            if arcpy.Exists(temp_table):
                arcpy.Delete_management(temp_table)
            arcpy.da.NumPyArrayToTable(struct[valid_indices], temp_table)
            # Create table view of the time dimension from netCDF file
            temp_mdtable = r"md_table"
            arcpy.MakeNetCDFTableView_md(dataset_name, 'time', temp_mdtable, 'time')
            # Get time field from table view
            desc_fields = arcpy.ListFields(temp_mdtable)
            time_field = desc_fields[[field.name for field in \
                desc_fields].index('time')]
            # Create a table in the output location
            arcpy.CreateTable_management(*os.path.split(output_filename))
            # Add fields from temporary table except for time field, which
            #   we add from the MakeNetCDFTableView
            desc_fields = arcpy.ListFields(temp_table)
            for index, field in enumerate(desc_fields):
                if field.type != 'OID':
                    if field.name == 'time':
                        arcpy.AddField_management(output_filename, time_field.name, \
                            time_field.type, time_field.precision, time_field.scale, \
                            time_field.length, time_field.aliasName, \
                            time_field.isNullable, time_field.required, \
                            time_field.domain)
                    else:
                        arcpy.AddField_management(output_filename, field.name, \
                            field.type, field.precision, field.scale, field.length, \
                            field.aliasName, field.isNullable, field.required,
                            field.domain)
            arcpy.DeleteField_management(output_filename,"Field1")
            # Create list of datetime objects from table view
            time_arc = []
            with arcpy.da.SearchCursor(temp_mdtable, 'time') as cursor:
                for row in cursor:
                    time_arc.append(row[0])
            del cursor
            # Get time variable from netCDF dataset
            time_var = dataset.variable('time')[:].tolist()
            # Loop through entries from temporary table
            with arcpy.da.InsertCursor(output_filename, '*') as ins_cursor:
                new_time_ind = ins_cursor.fields.index('time')
                with arcpy.da.SearchCursor(temp_table, '*') as cursor:
                    time_ind = cursor.fields.index('time')
                    for row in cursor:
                        new_row = list(row)
                        new_row[new_time_ind] = \
                            time_arc[time_var.index(row[time_ind])]
                        ins_cursor.insertRow(tuple(new_row))
            del ins_cursor
            del cursor
        else:
            # Write table
            if arcpy.env.overwriteOutput:
                arcpy.Delete_management(output_filename)
            arcpy.da.NumPyArrayToTable(struct[valid_indices], output_filename)

        return
dataset/ETHPy150Open Esri/solutions-geoprocessing-toolbox/suitability/toolboxes/scripts/MultidimensionSupplementalTools/MultidimensionSupplementalTools/Scripts/mds/tools/multidimensional_zonal_statistics_as_table.py/MultidimensionalZonalStatisticsAsTable.zonal_statistics_as_table_for_discrete/swapped operands of `in`