def zonal_statistics_as_table_for_grid(self, parameters, messages, dataset):
        """Performs zonal statistics as table assuming dataset contains
        regularly gridded variables."""

        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)
        # Lower left hand corner of netCDF variable spatial slice
        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)

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

        # 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 not dim_item in spat_list:
                dim_prod *= dim_shape[dim_index]
                dim_slices.append(xrange(dim_shape[dim_index]))
                dim_dims.append((str(dim_item), numpy.int))
            else:
                dim_slices.append(Ellipsis)

        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)

        # 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
        index = 0
        valid_indices = []
        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]

            # 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)):
                    # Set mask to netCDF NoData 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)

                    # Structure results into row
                    struct_data = [dataset.variable(dim_name[it_ind])[it] for \
                        it_ind, it in enumerate(list(varslice)) if \
                        it != Ellipsis]
                    struct_data.append(z)
                    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 type_parameter.valueAsText in \
                                self.statistics_multiple_values:
                            struct_dtype += [(it, numpy.float) for it in list(
                                result.dtype.names)]
                        else:
                            struct_dtype.append((str(
                                type_parameter.valueAsText), numpy.float))
                        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_grid/swapped operands of `in`