MeteoInfoLab脚本示例：计算不同区域平均值

风格而才

王老师你好，我有个问题问一下您。我有一个东亚数据，现在想把中国的数据mskout出来，求平均值，但是结果出来平均值是零，麻烦王老师帮忙看一下，谢谢！前四行的数据读取没问题。用的一下程序：
f=addfile(r'D:/MODIS/all.nc ')
v_cod = f['MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean']
tdata = v_cod[:,:,:]
us = shaperead('D:/MeteoInfo_1.4.7R2/MeteoInfo/map/china.shp')
mdata = tdata.maskout(us)
tave = mdata.ave()

MeteoInfo

风格而才 发表于 2018-3-30 16:37
王老师你好，我有个问题问一下您。我有一个东亚数据，现在想把中国的数据mskout出来，求平均值，但是结果出 ...

你在Console里输入 f 回车，看看文件信息。

风格而才

MeteoInfo 发表于 2018-3-30 16:44
你在Console里输入 f 回车，看看文件信息。

这是信息。
Dimensions: 6
        lat = 47;
        lon = 74;
        time = 168;
        latv = 2;
        lonv = 2;
        nv = 2;
X Dimension: Xmin = 70.5; Xmax = 143.5; Xsize = 74; Xdelta = 1.0
Y Dimension: Ymin = 8.5; Ymax = 54.5; Ysize = 47; Ydelta = 1.0
Global Attributes:
        : nco_openmp_thread_number = 1
        : Conventions = "CF-1.4"
        : start_time = "2003-01-01T00:00:00Z"
        : end_time = "2003-01-31T23:59:59Z"
        : temporal_resolution = "monthly"
        : NCO = "\"4.5.3\""
        : title = "MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean (69.6094E_7.7344N_144.1406E_54.8438N)"
        : history = "Fri Nov 24 02:25:22 2017: ncatted -a title,global,o,c,MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean (69.6094E_7.7344N_144.1406E_54.8438N) -a valid_range,,d,, -O -o ./subsetted.MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean.20030101.69E_7N_144E_54N.nc ./subsetted.MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean.20030101.69E_7N_144E_54N.nc"
Variations: 8
        double MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean(time,lat,lon);
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: _FillValue = -9999.0
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: long_name = "Combined Cloud Optical Thickness: Mean of Daily Mean"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: units = "1"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Level_2_Pixel_Values_Read_As = "Real"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Included_Level_2_Nighttime_Data = "False"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Quality_Assurance_Data_Set = "Quality_Assurance_1km"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Statistic_Type = "Simple"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: QA_Byte = 0S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: QA_Useful_Flag_Bit = 0S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: QA_Value_Start_Bit = 1S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: QA_Value_Num_Bits = 2S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Data_Set = "Quality_Assurance_1km"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Byte = 2S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Value_Start_Bit = 0S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Value_Num_Bits = 3S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Category_Values = 2S, 3S, 4S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Aggregation_Valid_Category_Values = 1S, 2S, 3S, 4S
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Derived_From_Level_3_Daily_Data_Set = "Cloud_Optical_Thickness_Combined_Mean"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Weighting = "Pixel_Weighted"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: Weighted_Parameter_Data_Set = "Cloud_Retrieval_Fraction_Combined_Pixel_Counts"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: standard_name = "Cloud Optical Thickness"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: quantity_type = "Cloud Properties"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: product_short_name = "MYD08_M3"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: product_version = "6"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: coordinates = "time lat lon"
                MYD08_M3_6_Cloud_Optical_Thickness_Combined_Mean_Mean: _ChunkSizes = 1, 47, 74
        int datamonth(time);
                datamonth: units = "hours since 1800-1-1 00:00:00"
                datamonth: long_name = "Time"
                datamonth: standard_name = "time"
                datamonth: axis = "T"
        float lat(lat);
                lat: units = "degrees_north"
                lat: standard_name = "latitude"
                lat: bounds = "lat_bnds"
                lat: _ChunkSizes = 47
        double lat_bnds(lat,latv);
                lat_bnds: units = "degrees_north"
        float lon(lon);
                lon: units = "degrees_east"
                lon: standard_name = "longitude"
                lon: bounds = "lon_bnds"
                lon: _ChunkSizes = 74
        double lon_bnds(lon,lonv);
                lon_bnds: units = "degrees_east"
        int time(time);
                time: units = "hours since 1800-1-1 00:00:00"
                time: long_name = "Time"
                time: standard_name = "time"
                time: axis = "T"
        int time_bnds(time,nv);
                time_bnds: units = "seconds since 1970-01-01 00:00:00"
                time_bnds: _ChunkSizes = 1, 2

MeteoInfo

风格而才 发表于 2018-3-30 16:52
这是信息。
Dimensions: 6
        lat = 47;

你看看tdata的平均值是不是0

风格而才

MeteoInfo 发表于 2018-3-30 19:59
你看看tdata的平均值是不是0

问题已经解决了，我没有固定时间维。然后maskout那也没写正确，正确写法mdata = year.maskout(us.shapes()[0])。谢谢王老师的回复。

cnn丸子宁

请问站点区域平均怎么作

MeteoInfo

cnn丸子宁 发表于 2019-5-7 20:19
请问站点区域平均怎么作

可以用 rmaskout 函数删除某个区域外的站点数据，然后再平均：

1. f = addfile_micaps(r'D:/Temp/micaps/19050914.000')
   
2. data = f['Temperature'][:]
   
3. lon = f['Longitude'][:]
   
4. lat = f['Latitude'][:]
   
5. lchina = shaperead('china')
   
6. data, lon, lat = rmaskout(data, lon, lat, lchina.shapes())

f = addfile_micaps(r'D:/Temp/micaps/19050914.000')<br /> data = f['Temperature'][:]<br /> lon = f['Longitude'][:]<br /> lat = f['Latitude'][:]<br /> lchina = shaperead('china')<br /> data, lon, lat = rmaskout(data, lon, lat, lchina.shapes())

CHch

fn = 'D:\CH\meteoinfo\emis_huabei4km_20161201.ncf'
f = addfile(fn)
data = f['NO2'][:,0,:,:]
ave = mean(data, axis=0)
tdata = ave[:,:]
China = shaperead('E:\data\map\country\province_2004.shp')
China.addfield('temp', 'float')
geoshow('China', facecolor='lightgray', edge=False)

i = 0
for rpoly in China.shapes():
    name = China.cellvalue('NAME2004', i)
    mdata = tdata.maskout(rpoly)
    tave = mdata.ave()
    tmin = mdata.min()
    tmax = mdata.max()
    print name + ', Ave: %.2f, Min: %.2f, Max: %.2f' %(tave, tmin, tmax)
    i += 1
求问大神，我按您的程序进行了修改，我用的是ncf格式数据，是4维数据，我选了第0层，然后对时间维进行了平均赋值给了tdata，但是在下面循环时没有数值，大神看看哪里不对呢

MeteoInfo

CHch 发表于 2019-11-22 09:52
fn = 'D:\CH\meteoinfo\emis_huabei4km_20161201.ncf'
f = addfile(fn)
data = f['NO2'][:,0,:,:]

可能是应为你的数据并非经纬度投影，而地图数据是经纬度投影，二者不匹配造成的。可以将地图数据投影为和数据一致，然后在 maskout 。

China = shaperead ...
China.project(f.proj)

sfz13851581045

同问，请教，全国空气质量的站点数据，怎么求各个省份区域内的平均？使用脚本实例的code，提示NDarray没有mask属性