nc全球数据中国区域的提取

vay

个人需要从全球数据中提取中国的，找遍全网，没一个能用的。有个人发布的什么matlab的视频，说的天花乱坠，反正我下载下来看完基本什么用没有。
说下思路，新建空的nc文件，通过中国的0.5格点经纬度，去全球数据中找对应的，然后摘到新建的nc文件里去。
我自己的数据不是单个的变量，有五个，所以新建的变量比较多。中国格点放到txt里了

代码：
import numpy as np
from netCDF4 import Dataset



dataset = Dataset('E:/0China/waterBodiesChina.nc','w',format='NETCDF4_CLASSIC') #创建一个nc文件
lat = dataset.createDimension('lat',360)
lon = dataset.createDimension('lon',720)
time = dataset.createDimension('time',None)#相关维度的创建
#创建变量
latitudes = dataset.createVariable('latitude',np.float32,('lat'))
longitudes = dataset.createVariable('longitude',np.float32,('lon'))
times = dataset.createVariable('time',np.float32,('time'))
fracWaterInp = dataset.createVariable('fracWaterInp',np.float32,('lat','lon','time'))
resMaxCapInp = dataset.createVariable('resMaxCapInp',np.float32,('lat','lon','time'))
resSfAreaInp = dataset.createVariable('resSfAreaInp',np.float32,('lat','lon','time'))
waterBodyIds = dataset.createVariable('waterBodyIds',np.float32,('lat','lon','time'))
waterBodyTyp = dataset.createVariable('waterBodyTyp',np.float32,('lat','lon','time'))
#变量参数设置
latitudes.units = 'degree_north'
longitudes.units = 'degree_east'
times.units = 'days since  2000-01-01 00:00:00'
times.calendar = 'gregorian'
fracWaterInp.units = '-'
resMaxCapInp.units = '10^6 m3'
resSfAreaInp.units = '10^6 m2'
waterBodyIds.units = '-'
waterBodyTyp.units = '-'

#创建经纬度矩阵
lats = np.arange(89.75,-89.75-0.5,-0.5)#创建矩阵
lons = np.arange(-179.75,179.75+0.5,0.5)
latitudes[:] = lats
longitudes[:] = lons

def loadData(fileName):
    inFile = open(fileName,'r')

    LON = []
    LAT = []

    for line in inFile:
        position = line.split()#每行数据按照空格分开成两部分
        LON.append(position[0])
        LAT.append(position[1])

    return LON,LAT

LON,LAT = loadData('E:/0China/China_0.5.txt')

source = Dataset('E:/0China/waterBodies30min.nc')


A = np.zeros((360,720),dtype=float)#创建0矩阵
A[A==-0] =np.nan#全部替换成空值

#from datetime import datetime
import arrow
from netCDF4 import date2num

#添加时间
for count in range(11):#原始数据是从1901-2010，我需要2000-2010，所以最后赋值的时候出现count+100，原始数据的第100层等于新数据第一层

    date = arrow.Arrow(2000,1,1).shift(years = count)
    times_data = date.naive
    times[count] = date2num(times_data, units=times.units, calendar = times.calendar)


    fracWaterInp[:,:,count] = A
    resMaxCapInp[:,:,count] = A
    resSfAreaInp[:,:,count] = A
    waterBodyIds[:,:,count] = A
    waterBodyTyp[:,:,count] = A


     for lo in range (3877):
          LONX = (float(LON[lo])+0.25)*2+359#经纬度转nc矩阵行列号，如果用xarray可以不用转，我没有用那个，用的笨法子
          LATX = (89.75-float(LAT[lo]))*2
          fracWaterInp[LATX,LONX,count] = source.variables['fracWaterInp'][count+100,LATX,LONX]
          resMaxCapInp[LATX,LONX,count] = source.variables['resMaxCapInp'][count+100,LATX,LONX]
          resSfAreaInp[LATX,LONX,count] = source.variables['resSfAreaInp'][count+100,LATX,LONX]
          waterBodyIds[LATX,LONX,count] = source.variables['waterBodyIds'][count+100,LATX,LONX]
          waterBodyTyp[LATX,LONX,count] = source.variables['waterBodyTyp'][count+100,LATX,LONX]




dataset.close()

werewolf

搞这么复杂干嘛，cdo一行解决了

vay

不会使，只是自己这么摸索出来的

DrLeviathan

werewolf 发表于 2019-6-13 13:54
搞这么复杂干嘛，cdo一行解决了

cdo应该怎么操作？

禾栗儿

好复杂

slchengh

请问你这个数据是全球水库水体数据集合吗？哪里可以下载？