Резервная сетка от центральных координат к внешним (то есть угловым) координатам

Есть ли доступный метод экстраполяции угловых положений сетки (синие точки) из центральных позиций сетки (красные точки)?

Сетка, с которой я работаю, не является прямоугольной, поэтому регулярная двулинейная интерполяция не выглядит так, как это было бы лучше всего; хотя, это просто так, что в моем заговоре мои данные используют pyplot.pcolormesh() , поэтому, возможно, это не так важно.

введите описание изображения здесь

Пример данных сетки :

 import numpy as np lons = np.array([[ 109.93299681, 109.08091365, 108.18301276, 107.23602539], [ 108.47911382, 107.60397996, 106.68325946, 105.71386119], [ 107.06790187, 106.17259769, 105.23214707, 104.2436463 ], [ 105.69908292, 104.78633156, 103.82905363, 102.82453812]]) lats = np.array([[ 83.6484245 , 83.81088466, 83.97177823, 84.13098916], [ 83.55459198, 83.71460466, 83.87294803, 84.02950188], [ 83.4569054 , 83.61444708, 83.77022192, 83.92410637], [ 83.35554612, 83.51060313, 83.6638013 , 83.81501464]]) 

2 Solutions collect form web for “Резервная сетка от центральных координат к внешним (то есть угловым) координатам”

Я не знаю ни одного надежного метода matplotlib для выполнения того, что вы просите, но у меня может быть другое решение. Мне часто приходится заполнять / экстраполировать в области сетки, где мне не хватает информации. Для этого я использую программу Fortran, которую я компилирую с помощью F2PY (который поставляется с numpy), который создает его в модуль python. Предполагая, что у вас есть компилятор Intel Fortran, вы скомпилируете его с помощью следующей команды: f2py --verbose --fcompiler=intelem -c -m extrapolate fill.f90 . Вы можете вызвать программу из python (см. Здесь полный пример):

  import extrapolate as ex undef=2.0e+35 tx=0.9*undef critx=0.01 cor=1.6 mxs=100 field = Zg field=np.where(abs(field) > 50 ,undef,field) field=ex.extrapolate.fill(int(1),int(grdROMS.xi_rho), int(1),int(grdROMS.eta_rho), float(tx), float(critx), float(cor), float(mxs), np.asarray(field, order='Fortran'), int(grdROMS.xi_rho), int(grdROMS.eta_rho)) 

Программа решает уравнение Лапласа с граничными условиями Неймана (dA / dn = 0) в RECTANGULAR-координатах с помощью итерационного метода для заполнения разумных значений в точках сетки, содержащих значения типа «undef». Это отлично работает для меня, и, может быть, вы сочтете это полезным. Программа доступна на моей учетной записи github.

Это не очень элегантный подход, который я использовал с помощью pyproj чтобы сначала вычислить расстояние и азимут между точками (с pyproj.Geod.inv , затем интерполировать / экстраполировать этот угол на необходимое расстояние (с pyproj.Geod.fwd ) к позиции пси.

код :

 def calc_psi_coords(lons, lats): ''' Calcuate psi points from centered grid points''' import numpy as np import pyproj # Create Geod object with WGS84 ellipsoid g = pyproj.Geod(ellps='WGS84') # Get grid field dimensions ydim, xdim = lons.shape # Create empty coord arrays lons_psi = np.zeros((ydim+1, xdim+1)) lats_psi = np.zeros((ydim+1, xdim+1)) # Calculate internal points #-------------------------- for j in range(ydim-1): for i in range(xdim-1): lon1 = lons[j,i] # top left point lat1 = lats[j,i] lon2 = lons[j+1,i+1] # bottom right point lat2 = lats[j+1,i+1] # Calc distance between points, find position at half of dist fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*0.5) # Assign to psi interior positions lons_psi[j+1,i+1] = lon_psi lats_psi[j+1,i+1] = lat_psi # Caclulate external points (not corners) #---------------------------------------- for j in range(ydim): # Left external points #~~~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[j+1,2] # left inside point lat1 = lats_psi[j+1,2] lon2 = lons_psi[j+1,1] # left outside point lat2 = lats_psi[j+1,1] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[j+1,0] = lon_psi lats_psi[j+1,0] = lat_psi # Right External points #~~~~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[j+1,-3] # right inside point lat1 = lats_psi[j+1,-3] lon2 = lons_psi[j+1,-2] # right outside point lat2 = lats_psi[j+1,-2] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[j+1,-1] = lon_psi lats_psi[j+1,-1] = lat_psi for i in range(xdim): # Top external points #~~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[2,i+1] # top inside point lat1 = lats_psi[2,i+1] lon2 = lons_psi[1,i+1] # top outside point lat2 = lats_psi[1,i+1] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[0,i+1] = lon_psi lats_psi[0,i+1] = lat_psi # Bottom external points #~~~~~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[-3,i+1] # bottom inside point lat1 = lats_psi[-3,i+1] lon2 = lons_psi[-2,i+1] # bottom outside point lat2 = lats_psi[-2,i+1] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[-1,i+1] = lon_psi lats_psi[-1,i+1] = lat_psi # Calculate corners: #------------------- # top left corner #~~~~~~~~~~~~~~~~ lon1 = lons_psi[2,2] # bottom right point lat1 = lats_psi[2,2] lon2 = lons_psi[1,1] # top left point lat2 = lats_psi[1,1] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[0,0] = lon_psi lats_psi[0,0] = lat_psi # top right corner #~~~~~~~~~~~~~~~~~ lon1 = lons_psi[2,-3] # bottom left point lat1 = lats_psi[2,-3] lon2 = lons_psi[1,-2] # top right point lat2 = lats_psi[1,-2] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[0,-1] = lon_psi lats_psi[0,-1] = lat_psi # bottom left corner #~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[-3,2] # top right point lat1 = lats_psi[-3,2] lon2 = lons_psi[-2,1] # bottom left point lat2 = lats_psi[-2,1] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[-1,0] = lon_psi lats_psi[-1,0] = lat_psi # bottom right corner #~~~~~~~~~~~~~~~~~~~~ lon1 = lons_psi[-3,-3] # top left point lat1 = lats_psi[-3,-3] lon2 = lons_psi[-2,-2] # bottom right point lat2 = lats_psi[-2,-2] # Calc dist between points, find position at dist*2 from pos1 fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2) lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.) lons_psi[-1,-1] = lon_psi lats_psi[-1,-1] = lat_psi return lons_psi, lats_psi 

Пример изображения (вокруг Дании / южной Швеции) :

введите описание изображения здесь

Python - лучший язык программирования в мире.