首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >Python:如何以不同的空间分辨率合并两个不同的netCdf填充?

Python:如何以不同的空间分辨率合并两个不同的netCdf填充?
EN

Stack Overflow用户
提问于 2022-03-29 15:25:11
回答 1查看 140关注 0票数 0

有可能合并两个不同空间分辨率的netCDF文件吗?

我有两个数据集。

第一个是ESA土地覆被数据集,空间解析为300m为netCDF。

第一个是居住在意大利的人口,其空间分辨率为100m,从WorldPop转换为geoTIFF,我将其转换为netCDF

这就是我要做的

代码语言:javascript
复制
## 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数据

代码语言:javascript
复制
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

我得到了意大利的数据

代码语言:javascript
复制
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数据集从100m300m的和,在esaItaly的相同像素中重新采样。

这就是我试过的

代码语言:javascript
复制
wp_inter = wp.interp(lat=esaItaly["lat"], lon=esaItaly["lon"])

但人口总数却要少得多。

代码语言:javascript
复制
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
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2022-03-29 19:21:28

做这件事有很多种方法,但最简单的方法可能是gdalbuildvrt

使用gdalbuildvrt --无论是从命令行还是从Python库--并构建一个VRT数据集。确保在结尾处列出了最高分辨率的文件--如果存在重叠,则最终的dataset获胜。

记住使用[-resolution {highest|lowest|average|user}]选项。

拥有复合数据集后,请使用gdal_translate -- CLI或Python --将其合并为单一单一数据集,格式为首选格式。

不要试图自己去实现它--它比看起来要复杂得多。

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/71664793

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档