Search code examples
pythonpython-iris

how to categorize lat/lon to find nearest city


I am playing with iris (really neat!) and I have a list of cities lat/lons I am interested to see average temperature over time. I have netcdf files with air temperatures covering entire country. I would like to tag data points in a cube with lat/lons closest to my cities so then I can easily get values I need just for these cities, or export data just for these cities somewhere.

I imagine I need to use add_categorised_coord somehow? https://scitools.org.uk/iris/docs/latest/iris/iris/coord_categorisation.html#iris.coord_categorisation.add_categorised_coord

I will appreciate an example. Thanks!


Solution

  • Assuming you have a gridded dataset of air temperature, a better solution would be to interpolate the data to given coordinate points, instead of "tagging" data points in a cube.

    This can be done by looping over cities and their coordinates and using cube.interpolate() method. See https://scitools.org.uk/iris/docs/latest/userguide/interpolation_and_regridding.html#cube-interpolation-and-regridding for examples.

    A more optimised solution would be to interpolate the data to all city points at once using the trajectory module:

    import iris
    import iris.analysis.trajectory as itraj
    import numpy as np
    
    # Create some dummy data
    nx = 10
    ny = 20
    
    lonc = iris.coords.DimCoord(
        np.linspace(-5, 10, nx), units="degrees", standard_name="longitude"
    )
    latc = iris.coords.DimCoord(
        np.linspace(45, 55, ny), units="degrees", standard_name="latitude"
    )
    cube = iris.cube.Cube(
        np.random.rand(ny, nx),
        dim_coords_and_dims=((latc, 0), (lonc, 1)),
        standard_name="x_wind",
        units="m s^-1",
        attributes=dict(title="dummy_data"),
    )
    
    # Create a Nx2 array of city coordinates
    city_coords = np.array([[50.7184, -3.5339], [48.8566, 2.3522], [52.6401898, 1.2517445]])
    
    # Interpolate the data to the given points
    sample_points = [("latitude", city_coords[:, 0]), ("longitude", city_coords[:, 1])]
    cube_values_in_cities = itraj.interpolate(cube, sample_points, "linear")
    

    Hope this helps.