<jupyter_start><jupyter_text>Read and visualise data<jupyter_code>date = iris.time.PartialDateTime(year=2006, month=6, day=22)
time_constraint = iris.Constraint(time=date)
level_constraint = iris.Constraint(lev=5)

with iris.FUTURE.context(cell_datetime_objects=True):
    u_cube = iris.load_cube(u_file, 'eastward_wind' & time_constraint)  
    v_cube = iris.load_cube(v_file, 'northward_wind' & time_constraint) 
    #thetao_cube = iris.load_cube(thetao_file, 'sea_water_potential_temperature')
%matplotlib inline

def plot_winds(u, v):
    """Plot the wind."""

    proj = ccrs.PlateCarree(central_longitude=-180.0)
    fig = plt.figure(figsize=(12, 5))

    ax = plt.subplot(121, projection=proj)
    qplt.contourf(u, levels=ticks, cmap=cmap)
    plt.gca().coastlines()

    ax = plt.subplot(122, projection=proj)
    qplt.contourf(v, levels=ticks, cmap=cmap)
    plt.gca().coastlines()

    iplt.show()
plot_winds(u_cube, v_cube)<jupyter_output><empty_output><jupyter_text>Rotate the wind<jupyter_code>rotated_cs = iris.coord_systems.RotatedGeogCS(20, 260)
urot_cube, vrot_cube = iris.analysis.cartography.rotate_winds(u_cube, v_cube, rotated_cs)
plot_winds(urot_cube, vrot_cube)
print vrot_cube
vrot_cube.coords()
x_values = vrot_cube.coord('projection_x_coordinate').points
y_values = vrot_cube.coord('projection_y_coordinate').points
lats = vrot_cube.coord('latitude').points
lons = vrot_cube.coord('longitude').points

print x_values.min(), x_values.max()
print y_values.min(), y_values.max()
x_values_360 = numpy.where(x_values < 0, x_values + 360, x_values)
print x_values_360.min(), x_values_360.max()
print lats[10]
print lons[65]
print 'new lon:', x_values[10, 65]
print 'new lat:', y_values[10, 65]
lon_mesh, lat_mesh = numpy.meshgrid(lons, lats)

lat_flat = lat_mesh.flatten()
lon_flat = lon_mesh.flatten()

x_flat = x_values.flatten()
x_flat_360 = x_values_360.flatten()
y_flat = y_values.flatten()
len(lat_flat)
len(x_flat)
def convert_point_1d(orig_lat, orig_lon):
    
    lat_indexes = numpy.where(lat_flat == orig_lat)
    lon_indexes = numpy.where(lon_flat == orig_lon)

    index = numpy.intersect1d(lat_indexes[0], lon_indexes[0])
    assert len(index) == 1
    
    print 'rotated lat (1D method):', y_flat[index]
    print 'rotated lon (1D method):', x_flat[index]
    
def convert_point_2d(orig_lat, orig_lon):
    
    lat_index = numpy.where(lats == orig_lat)
    lon_index = numpy.where(lons == orig_lon)
    
    print 'rotated lat (2D method):', y_values[lat_index, lon_index]
    print 'rotated lon (2D method):', x_values[lat_index, lon_index]
test_lat = 82.5
test_lon = 130.5

convert_point_1d(test_lat, test_lon)
convert_point_2d(test_lat, test_lon)<jupyter_output>rotated lat (1D method): [ 15.13393305]
rotated lon (1D method): [ 5.98889898]
rotated lat (2D method): [[ 15.13393305]]
rotated lon (2D method): [[ 5.98889898]]<jupyter_text>Regrid the data, scipy methodMy search path goes from 115E to 225E (or 135W), which means it straddles 180E. To avoid any edge issues, I'm going to want to do the regridding with longitude set from 0 to 360.<jupyter_code>points = numpy.column_stack((x_flat_360, y_flat))
values = vrot_cube.data.flatten()

grid_lons, grid_lats = numpy.meshgrid(numpy.arange(0, 360), numpy.arange(-90, 91))                           


regridded_data = scipy.interpolate.griddata(points, values, (grid_lons, grid_lats), 
                                            method='linear', fill_value=0)
print regridded_data.shape
print values.max(), values.min()
print regridded_data.max(), regridded_data.min()
regridded_data_clean = numpy.where(regridded_data < values.min(), values.min(), regridded_data) 
regridded_data_clean = numpy.where(regridded_data_clean > values.max(), values.max(), regridded_data_clean)

import matplotlib.pyplot
image = matplotlib.pyplot.imshow(regridded_data_clean)
matplotlib.pyplot.show(image)<jupyter_output><empty_output><jupyter_text>Not sure why this is upsidedown?<jupyter_code>regridded_data_clean.shape<jupyter_output><empty_output><jupyter_text>Creating a cube from scratchThe documentation is [here](http://scitools.org.uk/iris/docs/latest/iris/iris/cube.html)<jupyter_code>from iris.coords import DimCoord
from iris.cube import Cube

latitude = DimCoord(numpy.arange(-90, 91, 1),
                    standard_name='latitude',
                    units='degrees')
longitude = DimCoord(numpy.arange(0, 360, 1),                     
                     standard_name='longitude',
                     units='degrees')
new_cube = Cube(regridded_data_clean,
                dim_coords_and_dims=[(latitude, 0),
                                     (longitude, 1)])

new_cube.long_name = vrot_cube.long_name
new_cube.units = vrot_cube.units
print new_cube<jupyter_output>transformed_northward_wind / (m s**-1) (latitude: 181; longitude: 360)
     Dimension coordinates:
          latitude                              x               -
          longitude                             -               x<jupyter_text>As discussed in [this issue](https://github.com/SciTools/cartopy/issues/644) the specified projection for the input data should not specify a `central_rotated_longitude`, even if it is 0 to 360 rather than -180 to 180...<jupyter_code>out_proj = ccrs.RotatedPole(260, 20, central_rotated_longitude=180)
fig = plt.figure(figsize=(12, 5))

ax = plt.subplot(111, projection=out_proj)
x = new_cube.coord('longitude').points
y = new_cube.coord('latitude').points

ticks = numpy.arange(-32, 36, 4)
cf = ax.contourf(x, y, new_cube.data, levels=ticks, cmap=cmap,
                 transform=ccrs.RotatedPole(260, 20))
plt.colorbar(cf)
plt.gca().coastlines()

iplt.show()<jupyter_output>WARNING:shapely.geos:Self-intersection at or near point 179.99999999999997 79.000000362338554<jupyter_text>Regrid the data, iris method<jupyter_code>vrot_cube.coords('latitude')[0].coord_system = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
vrot_cube.coords('longitude')[0].coord_system = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
def make_grid():
    """Make an empty Cube."""
    
    coordsys = iris.coord_systems.RotatedGeogCS(20, 260)
    lat_values = numpy.arange(-90, 91, 1)
    lon_values = numpy.arange(0, 360, 1)
    
    latitude = DimCoord(lat_values,
                        standard_name='latitude',
                        units='degrees_north',
                        coord_system=coordsys)
    longitude = DimCoord(lon_values,                     
                         standard_name='longitude',
                         units='degrees_east',
                         coord_system=coordsys)

    dummy_data = numpy.zeros((len(lat_values), len(lon_values)))
      
    new_cube = Cube(dummy_data, dim_coords_and_dims=[(latitude, 0), (longitude, 1)])
    
    return new_cube
target_grid_cube = make_grid()
scheme = iris.analysis.Linear()  # Or Nearest()
vrot_iris_regrid = vrot_cube.regrid(target_grid_cube, scheme)
print vrot_iris_regrid
out_proj = ccrs.RotatedPole(260, 20, central_rotated_longitude=180)
fig = plt.figure(figsize=(12, 5))

ax = plt.subplot(111, projection=out_proj)
x = vrot_iris_regrid.coord('longitude').points
y = vrot_iris_regrid.coord('latitude').points

ticks = numpy.arange(-32, 36, 4)
cf = ax.contourf(x, y, vrot_iris_regrid.data, levels=ticks, cmap=cmap,
                 transform=ccrs.RotatedPole(260, 20))
plt.colorbar(cf)
plt.gca().coastlines()

iplt.show()<jupyter_output>WARNING:shapely.geos:Self-intersection at or near point 179.99999999999997 78.220881880335725
{'fraction_non_alphanumeric': 0.12013805920615957, 'fraction_numerical': 0.041550511084561266, 'mean_word_length': 4.012641383898869, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 26, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '434670', 'license': 'MIT', 'n_tokens_mistral': 3049, 'n_tokens_neox': 2769, 'n_words': 484, 'path': 'development/wind_rotation.ipynb', 'repo_name': 'DamienIrving/climate-analysis', 'size': '7533'}