Partial Dependence Plot orients differently compared to Partial Dependence values

### Describe the bug

The issue is that the 2D partial dependence plot from scikit-learn orients in a different way that what you would get using raw pdp values from sklearn as well.



### Steps/Code to Reproduce

```python
import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingRegressor


DATA_COUNT = 10000

RANDOM_SEED = 137

TRAIN_FRACTION = 0.9

np.random.seed(RANDOM_SEED)

xs = pd.DataFrame(np.random.uniform(size = (DATA_COUNT, 3)))
xs.columns = ['x0', 'x1', 'x2']

y = pd.DataFrame(xs.x0*xs.x1 + xs.x2/xs.x1 + pd.Series(0.1*np.random.randn(DATA_COUNT)))
y.columns = ['y']

gbr_1 = GradientBoostingRegressor(random_state = RANDOM_SEED)
gbr_1.fit(xs, y)


from sklearn.inspection import partial_dependence, PartialDependenceDisplay

PartialDependenceDisplay.from_estimator(gbr_1, xs, [(0, 1)],grid_resolution=15,percentiles=(0.01,0.99),
                                            contour_kw=dict(cmap='bwr'))

# Versus 

pdp1 = partial_dependence(gbr_1, xs, features=[(0, 1)], kind='average', percentiles=(0.01,0.99),grid_resolution=15)
y_pred1 = pdp1['average']
grid1 = pdp1['values']

x,y = np.meshgrid(grid1[0],grid1[1])
plt.contourf(x,y,y_pred1[0],cmap='bwr',levels=7)


```

### Expected Results

```python
from sklearn.inspection import plot_partial_dependence,partial_dependence, PartialDependenceDisplay

PartialDependenceDisplay.from_estimator(gbr_1, xs, [(0, 2)],grid_resolution=15,percentiles=(0.01,0.99),
                                            contour_kw=dict(cmap='bwr'))
```


![image](https://user-images.githubusercontent.com/49596459/214959496-168287b8-6be1-447f-ae39-3b536c9ace4e.png)


### Actual Results

```python
pdp1 = partial_dependence(gbr_1, xs, features=[(0, 2)], kind='average', percentiles=(0.01,0.99),grid_resolution=15)
y_pred1 = pdp1['average']
grid1 = pdp1['values']

x,y = np.meshgrid(grid1[0],grid1[1])
plt.contourf(x,y,y_pred1[0],cmap='bwr',levels=7)
plt.xlabel('x0')
plt.ylabel('x2')
```

![image](https://user-images.githubusercontent.com/49596459/214959548-c40d8324-c204-40d4-ad7f-0d2cee1e1f70.png)


### Versions

```shell
System:
    python: 3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 07:04:59) [GCC 10.3.0]
executable: /glade/work/chayan/conda-envs/pangeo_chayan/bin/python
   machine: Linux-3.10.0-1127.18.2.el7.x86_64-x86_64-with-glibc2.17

Python dependencies:
      sklearn: 1.1.1
          pip: 22.1.2
   setuptools: 63.2.0
        numpy: 1.23.5
        scipy: 1.8.1
       Cython: 0.29.30
       pandas: 1.4.3
   matplotlib: 3.5.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /glade/work/chayan/conda-envs/pangeo_chayan/lib/libopenblasp-r0.3.20.so
        version: 0.3.20
threading_layer: pthreads
   architecture: SkylakeX
    num_threads: 1

       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /glade/work/chayan/conda-envs/pangeo_chayan/lib/libgomp.so.1.0.0
        version: None
    num_threads: 1

       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /glade/work/chayan/conda-envs/pangeo_chayan/lib/python3.10/site-packages/xgboost.libs/libgomp-a34b3233.so.1.0.0
        version: None
    num_threads: 1
```

COMMENT:
The issue is how you are using `contourf`: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.contourf.html

From the documentation, it states that `Z` (here the partial dependence values) should be of shape (M, N) where M is `len(Y)` (second feature) and N is `len(X)` (first feature). Then, you need to transpose.

If you want to convince yourself, you can have a look at a simple example, with non-equal shape for X and Y:

```python
import matplotlib.pyplot as plt
import numpy as np

In [46]: xx, yy = np.meshgrid([0, 1, 2], [3, 4])

In [47]: z = np.arange(0, 60, 10).reshape(3, 2)

In [48]: xx
Out[48]: 
array([[0, 1, 2],
       [0, 1, 2]])

In [49]: yy
Out[49]: 
array([[3, 3, 3],
       [4, 4, 4]])

In [50]: z
Out[50]: 
array([[ 0, 10],
       [20, 30],
       [40, 50]])

In [51]: plt.contourf(xx, yy, z.T)
Out[51]: <matplotlib.contour.QuadContourSet at 0x1278098d0>

In [52]: plt.colorbar()
Out[52]: <matplotlib.colorbar.Colorbar at 0x1277b7400>
```

`z` here if of shape `N, M` as we store the `pdp_values`. If we don't transpose then we cannot plot.

Labels: Bug, Needs Triage