How to apply a point transformation to many points?

user_10732424 0
Viewed 97

I have a gridded temperature dataset and a list of weather stations across the country and their latitudes and longitudes. I want to find the grid points that are nearest to the weather stations. My gridded data has coordinates x,y which latitude and longitude are a function of. enter image description here

I found that the simplest way of finding the nearest grid point is to first transform the latitude and longitude (Lat, Lon) of the weather stations to x and y values and then find the nearest grid point. I did that for one station (lat= , lon= ) by doing the following:

import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset
import numpy as np
from cartopy import config
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import xarray as xr
import pandas as pd
import netCDF4 as nc

#open gridded data
df=xr.open_dataset('/home/mmartin/LauNath/air.2m.2015.nc')

#open weather station data
CMStations=pd.read_csv('Slope95.csv')

import cartopy.crs as ccrs

# Example - your x and y coordinates are in a Lambert Conformal projection
data_crs = ccrs.LambertConformal(central_longitude=-107.0,central_latitude=50.0,standard_parallels = (50, 50.000001),false_easting=5632642.22547,false_northing=4612545.65137)

# Transform the point - src_crs is always Plate Carree for lat/lon grid
x, y = data_crs.transform_point(-94.5786,39.0997, src_crs=ccrs.PlateCarree())

# Now you can select data
ks=df.sel(x=x, y=y, method='nearest')

How would I apply this to all of the weather stations latitudes and longitudes (Lat,Lon)?

cartopy python-xarray pandas numpy python
2 Answers

You can create a geopandas GeoDataFrame from x, y columns using geopandas.points_from_xy. I'll assume these points are WGS84/EPSG4326:

import geopandas as gpd

stations = gpd.GeoDataFrame(
    CMStations,
    geometry=gpd.points_from_xy(
        CMStations.Lon, CMStations.Lat, crs="epsg:4326"  # assume WGS84
    ),
)

Now, we can use geopandas.GeoDataFrame.to_crs to transform all the points at once:

stations_xy = stations.to_crs(data_crs)

Finally, we can use xarray's advanced indexing, using DataArrays with lat/lon data and station ID as coordinates, to reshape the x/y data to the shape of the CMStations index:

station_x = stations.geometry.x.to_xarray()
station_y = stations.geometry.y.to_xarray()

# use these to select from xarray Dataset ds
station_data = ds.sel(y=station_y, x=station_x, method="nearest")

If desired, you could set a station ID column to be the index first with CMStations.set_index("station_id") to get the station_id column as the dataset dimension which replaces x and y.

user_38887190

There is no need to use geopandas in here... just use crs.transform_points() instead of crs.transform_point() and pass the coordinates as arrays!

import numpy as np
import cartopy.crs as ccrs

data_crs = ccrs.LambertConformal(central_longitude=-107.0,central_latitude=50.0,standard_parallels = (50, 50.000001),false_easting=5632642.22547,false_northing=4612545.65137)

lon, lat = np.array([1,2,3]), np.array([1,2,3])
data_crs.transform_points(ccrs.PlateCarree(), lon, lat)

which will return an array of the projected coordinates:

array([[16972983.1673108 ,  8528848.37931063,        0.        ],
       [16841398.80456616,  8697676.02704447,        0.        ],
       [16709244.32834945,  8862533.81411212,        0.        ]])

... also... if you really have a lot of points to transform (and maybe use some crs not yet supported by cartopy) you might want to have a look at PyProj directly since it provides a lot more functionality and also some tricks to speed up transformations. (it's used under the hood by cartopy as well so you should already have it installed!)

user_97034510
Related