Python:如何合并两个具有不同空间分辨率的不同netCdf填充?
是否可以将两个NetCDF文件合并为不同的空间分辨率?
我有两个数据集。
第一个是 esa land Cover cover dataset 300m
作为NETCDF的重置。
第一个是居住在意大利的人口,其空间分辨率为100m
来自 worldpop 作为geotiff
,我将其转换为netcdf
。
这是我正在这样做的
## convert GeoTiff to netCDF
from osgeo import gdal
fin = ita_ppp_2000.tif'
fout = 'ita_ppp_2000.nc'
ds = gdal.Translate(fout, fin, format='NetCDF')
我下载了ESA CCI数据
year = 2000
import cdsapi
c.retrieve(
'satellite-land-cover',
{
'variable': 'all',
'format': 'zip',
'version': 'v2.1.1',
'year': str(y),
},
'download_%d.zip'%y) ## we need to unzip it
fn = 'ESACCI-LC-L4-LCCS-Map-300m-P1Y-%d-v2.0.7cds.nc'%year## Global 300m resolution
,我收到了Italy New的数据,
def clipEsa(fn,x0,x1,y0,y1):
dnc = xr.open_dataset(fn)
lat = dnc.variables['lat'][:]
lon = dnc.variables['lon'][:]
# All indices in bounding box:
where_j = np.where((lon >= x0) & (lon <= x1))[0]
where_i = np.where((lat >= y0) & (lat <= y1))[0]
# Start and end+1 indices in each dimension:
i0 = where_i[0]
i1 = where_i[-1]+1
j0 = where_j[0]
j1 = where_j[-1]+1
longitude = dnc["lccs_class"]["lon"].values[j0:j1]
latitude = dnc["lccs_class"]["lat"].values[i0:i1]
time = dnc['lccs_class']['time'][0]
return dnc.sel(time=time, lon=longitude, lat=latitude)
wp = xr.open_dataset(fout) ## Italian population with 100m resolution
bb = [wp.lon[0].values, wp.lon[-1].values, wp.lat[0].values, wp.lat[-1].values] ## bounding box
esaItaly = clipEsa(fn,bb[0],bb[1],bb[2],bb[3]) ## ESA CCI clipped for Italy
我希望将两个数据集都以300m
的空间分辨率为单位。特别是我想用一个和将wp
数据集从100m
到300m
中的数据集中,以sysaitaly >
这是我尝试的
wp_inter = wp.interp(lat=esaItaly["lat"], lon=esaItaly["lon"])
,但人口总数要低得多。
sum(wp_inter['Band1'].values[wp_inter['Band1'].values>0])
5038174.5 ## population interpolated
sum(wp.Band1.values[wp.Band1.values>0])
56780870.0 ## original population
Is it possible to merge two netCDF files with different spatial resolution?
I have two datasets.
The first one is the ESA Land Cover dataset with a spatial resoltion of 300m
as netCDF.
The first one is the population living in Italy with a spatial resolution of 100m
from WorldPop as a geoTIFF
that I convert as netCDF
.
This is what I am doing
## convert GeoTiff to netCDF
from osgeo import gdal
fin = ita_ppp_2000.tif'
fout = 'ita_ppp_2000.nc'
ds = gdal.Translate(fout, fin, format='NetCDF')
I download the ESA CCI data
year = 2000
import cdsapi
c.retrieve(
'satellite-land-cover',
{
'variable': 'all',
'format': 'zip',
'version': 'v2.1.1',
'year': str(y),
},
'download_%d.zip'%y) ## we need to unzip it
fn = 'ESACCI-LC-L4-LCCS-Map-300m-P1Y-%d-v2.0.7cds.nc'%year## Global 300m resolution
I get the data for Italy
def clipEsa(fn,x0,x1,y0,y1):
dnc = xr.open_dataset(fn)
lat = dnc.variables['lat'][:]
lon = dnc.variables['lon'][:]
# All indices in bounding box:
where_j = np.where((lon >= x0) & (lon <= x1))[0]
where_i = np.where((lat >= y0) & (lat <= y1))[0]
# Start and end+1 indices in each dimension:
i0 = where_i[0]
i1 = where_i[-1]+1
j0 = where_j[0]
j1 = where_j[-1]+1
longitude = dnc["lccs_class"]["lon"].values[j0:j1]
latitude = dnc["lccs_class"]["lat"].values[i0:i1]
time = dnc['lccs_class']['time'][0]
return dnc.sel(time=time, lon=longitude, lat=latitude)
wp = xr.open_dataset(fout) ## Italian population with 100m resolution
bb = [wp.lon[0].values, wp.lon[-1].values, wp.lat[0].values, wp.lat[-1].values] ## bounding box
esaItaly = clipEsa(fn,bb[0],bb[1],bb[2],bb[3]) ## ESA CCI clipped for Italy
New I would like to have both the datasets at the spatial resolution of 300m
. In particular I would like to resample with a sum the wp
dataset from 100m
to 300m
in the same pixels of esaItaly
This is what I tried
wp_inter = wp.interp(lat=esaItaly["lat"], lon=esaItaly["lon"])
but the total amount of population is much lower.
sum(wp_inter['Band1'].values[wp_inter['Band1'].values>0])
5038174.5 ## population interpolated
sum(wp.Band1.values[wp.Band1.values>0])
56780870.0 ## original population
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(1)
有很多方法可以做到这一点,但最简单的方法是
gdalbuildvrt
。使用
gdalbuildvrt
- 从python库中的命令行 - 都可以构建vrt
dataset。确保最高的分辨率文件在结束时列出 - 如果最终数据集胜利重叠。请记住使用
[ - 分辨率{最高|最低|平均| user}]
选项。拥有复合数据集后,请使用
gdal_translate
-cli或python-将其以您首选的格式合并到单个单片数据集中。不要尝试自己实施 - 它比看起来更复杂。
There are many ways to do it, but probably the easiest one is by
gdalbuildvrt
.Use
gdalbuildvrt
- either from the command-line either from the Python library - and build aVRT
dataset. Make sure the highest resolution files are listed towards the end - if there is overlapping the final dataset wins.Remember to use
[-resolution {highest|lowest|average|user}]
option.Once you have a composite Dataset, use
gdal_translate
- CLI or Python - to consolidate it to a single monolithic Dataset in your preferred format.Don't try to implement this yourself - it is more complicated than it might seem.