Рассчитать среднее значение переменных в селективной области, в файле netCDF с сеткой
Допустим, у нас есть данные об осадках TRMM, каждый файл представляет данные за каждый месяц. Например, файлы в папке:
3B42.1998.01.01.7A.nc,
3B42.1998.02.01.7A.nc,
3B42.1998.03.01.7A.nc,
3B42.1998.04.01.7A.nc,
3B42.1998.05.01.7A.nc,
......
......
3B42.2010.11.01.7A.nc,
3B42.2010.12.01.7A.nc.
Эти файлы имеют следующие размеры: Xsize=1440, Ysize=400, Zsize=1,Tsize=1. Долгота установлена от 0 до 360, а широта установлена от -50 до 50. Я хочу рассчитать количество осадков в определенном регионе, скажем, между lon=98.5, lon=100 and lat=4, lat=6.5, Это означает, что читать переменные только в этом регионе -:
--------------------
|lon:98.5 lat:6.5|
| |
|lat:4 lon:100 |
---------------------
Я делал это в GrADS (система анализа и отображения сетки). В GrADS это можно сделать так: (упрощенная версия)
yy=1998
while yr < 2011
'sdfopen f:\data\trmm\3B42.'yy'.12.01.7A.nc'
'd aave(pcp,lon=98.5,lon=100.0,lat=4.0,lat=6.5)'
res=subwrd(result,4)
rec=write('d:\precip.sp.TRMM3B42.1.'yy'.csv',res,append)
yy = yy+1
endwhile
Я пытался сделать то же самое в Python, но что-то пошло не так. После нескольких предложений вот и я:
import csv
import netCDF4 as nc
import numpy as np
#calculating december only
f = nc.MFDataset('d:/data/trmm/3B43.????.12.01.7A.nc')#maybe I shouldn't do MFDataset?
pcpt = f.variables['pcp']
lon = f.variables['longitude']
lat = f.variables['latitude']
# Determine which longitudes
latidx1 = (lat >=4.0 ) & (lat <=6.5 )
lonidx1 = (lon >=98.5 ) & (lon <=100.0 )
rainf1 = pcpt[:]
rainf1 = rainf1[:, latidx1][..., lonidx1]
rainf_1 = rainf1
with open('d:/trmmtest.csv', 'wb') as fp:
a = csv.writer(fp)
for i in rainf_1:
a.writerow([i])
Этот скрипт создает список (в моем случае) 15 значений в CSV-файле. Но когда я пытаюсь получить значения для другого региона и настроить то, что считаю необходимым, скажем:
latidx2 = (lat >=1.0 ) & (lat <=1.5 )
lonidx2 = (lon >=102.75 ) & (lon <=103.25 )
rainf2 = pcpt[:]
rainf2 = rainf2[:, latidx2][..., lonidx2]
rainf_2 = rainf2
Я получаю те же значения, что и первый.
firstarea=[0.511935,1.0771,0.613548,1.48839,0.445161,1.39161,1.03548,0.452903, 3.07725,2.84613 0.701613,2.10581,2.47839,3.84097,2.41065,1.38387]
secondarea=[0.511935,1.0771,0.613548,1.48839,0.445161,1.39161,1.03548,0.452903, 3.07725,2.84613,0.701613,2.10581,2.47839,3.84097,2.41065,1.38387]
Я проводил тестирование по отдельным сценариям, они по-прежнему дают те же значения Я проверил на карте (построенной ранее), значения в этих двух регионах разные (в среднем за декабрь).
Есть идеи почему? Есть ли другой элегантный способ написать это? Спасибо.
4 ответа
Через некоторое время мне снова удалось взглянуть на эту проблему, и, по-видимому, описанный выше метод почти верен. После нескольких настроек, протестированных в одном файле данных и перепроверенных с помощью решения GrADS, я получил что-то вроде этого:
f = nc.Dataset('~/data/TRMM3H/3B42.19980101.12.7A.nc')
pcpt = f.variables['pcp'][:]
lon = f.variables['longitude'][:]
lat = f.variables['latitude'][:]
#select two regions
latidx1 = (lat >=4. ) & (lat <=6.5 )
lonidx1 = (lon >=100.5 ) & (lon <=101.5 )
latidx2 = (lat >=2.5 ) & (lat <=5.0 )
lonidx2 = (lon >=101. ) & (lon <=102. )
rainf = pcpt[:]
#these basically listing the values in an array (2 in this case)
rainf1 = rainf[:, latidx1][..., lonidx1]
rainf2 = rainf[:, latidx2][..., lonidx2]
rainf_1 = rainf1
rainf_2 = rainf2
#time to get the mean values
print np.mean(rainf_1)
print "............."
print np.mean(rainf_2)
print "............."
это дало мне такие результаты:
>>> execfile('find_percentile.py')
0.7830327034
.............
1.56235361099
.............
Результаты совпадают при расчете с GrADS.
Отредактировано после предложения:
f = nc.Dataset('~/data/TRMM3H/3B42.19980101.12.7A.nc')
pcpt = f.variables['pcp'][:]
lon = f.variables['longitude'][:]
lat = f.variables['latitude'][:]
#select two regions
latidx1 = (lat >=4. ) & (lat <=6.5 )
lonidx1 = (lon >=100.5 ) & (lon <=101.5 )
latidx2 = (lat >=2.5 ) & (lat <=5.0 )
lonidx2 = (lon >=101. ) & (lon <=102. )
#these basically listing the values in an array (2 in this case)
rainf1 = pcpt[:, latidx1][..., lonidx1]
rainf2 = pcpt[:, latidx2][..., lonidx2]
rainf_1 = rainf1
rainf_2 = rainf2
#time to get the mean values
print np.mean(rainf_1)
print "............."
print np.mean(rainf_2)
print "............."
Вернемся к исходному вопросу, делая это в нескольких файлах и печатая его в файле txt/csv, все еще в стадии разработки (и тестирования).
Я просто хочу отметить, что решение Fir Nor неверно, так как вы не можете просто использовать среднее арифметическое (np.mean) при работе с пространственными данными на регулярной решетке широта / долгота, как здесь, так как сетка размер клетки меняется по мере движения к полюсам!
Куда лучше не беспокоиться об этом и делать операцию с CDO:
cdo fldmean -sellonlatbox,98.5,100,4.5,6 3B42.1998.05.01.7A.nc boxav.nc
Если вы работаете в Linux, эту проблему можно решить с помощью nctoolkit (http://nctoolkit.readthedocs.io/en/latest/). Все должно делать следующее:
import nctoolkit as nc
ff = '~/data/TRMM3H/3B42.19980101.12.7A.nc'
data = nc.open_data(ff)
data.clip(lon = [98.5, 100], lat = [4, 6.5])
data.spatial_mean()
Примечание: здесь в качестве бэкэнда используется CDO, а Space_mean вычисляет среднее значение, взвешенное по площади каждой ячейки сетки.
Я считаю, что это легко сделать с помощью пакета easymore.
Первый шаг - это создание шейп-файла. Это может быть любая форма (например, точка, подбассейн или прямоугольник). В вашем случае это будет прямоугольный шейп-файл с одной формой, определяющей границу. Это можно сделать в QGIS, ArcGIS или на Python:
Создание файла формы из списка координат ограничивающей рамки
Затем нужно вызвать пакет easymore python и переназначить переменную в интересующий шейп-файл, выполнив следующие простые действия:
# loading EASYMORE
from easymore.easymore import easymore
# initializing EASYMORE object
esmr = easymore()
# specifying EASYMORE objects
# name of the case
esmr.case_name = 'CMIP6_3B43'
# temporary path that the EASYMORE generated GIS files and remapped file will be saved
esmr.temp_dir = 'path/temporary/'
# name of target shapefile that the source netcdf files should be remapped to;
# it was created in the first step
esmr.target_shp = 'path/target_shapefiles/box.shp'
# name of netCDF file(s); multiple files can be specified with *
esmr.source_nc = ' d:/data/trmm/3B43*.nc'
# name of variables from source netCDF file(s) to be remapped
esmr.var_names = ['pcp']
# name of variable longitude in source netCDF files
esmr.var_lon = 'longitude'
# name of variable latitude in source netCDF files
esmr.var_lat = 'latitude'
# name of variable time in source netCDF file; should be always time
esmr.var_time = 'time'
# location where the remapped netCDF, csv file will be saved
esmr.output_dir = 'path/output/'
# if required that the remapped values to be saved as csv as well
esmr.save_csv = True
# execute EASYMORE nc remapper
esmr.nc_remapper()
Этот код будет генерировать как переназначенные файлы nc, так и его версию csv в выходном каталоге для каждого исходного файла nc. Переназначенные файлы будут представлять собой среднее по площади количество осадков по интересующей форме (ам) с исходным разрешением по времени (например, день). Затем вы можете легко увеличить их до ежемесячного временного шага и сравнить.
Преимущества:
1- С помощью этого пакета вы можете предоставить шейп-файл с несколькими формами (интересующими областями), который выполняет переназначение за один раз. В качестве примера вы можете вместо этого предоставить шейп-файл стран мира.
2- Если ваше поле меньше, чем сетка (многоугольник) или точка, возвращаемое значение будет сеткой, в которой находится маленький прямоугольник или точка.
3- Переназначение и взвешивание выполняются в равной области для учета различных сеток одинаковой площади в WGS84 в более высоких латах.
4- Код достаточно умен, поэтому вам не нужно беспокоиться о формате от 0 до 360 lon и о целевом шейп-файле в формате от -180 до 180 lon. Например, если поле находится в NA, отрицательное значение lon, шейп-файл может быть задан в отрицательном формате lon от -180 до 180, а файл NC - с неотрицательными значениями lon (от 0 до 360).
Страница GitHub с дополнительными примерами: