MeteoInfoLab脚本示例：Sea ice netCDF数据

real-slimshady

MeteoInfo 发表于 2016-4-1 20:41
原始数据就是这样呀，没有特殊处理。图例用cmap参数控制。

老师，如何让图像顺时针或逆时针旋转90°？还有我用了add point但是无法添加点

MeteoInfo

real-slimshady 发表于 2016-4-5 16:08
老师，如何让图像顺时针或逆时针旋转90°？还有我用了add point但是无法添加点

顺指针旋转90度，实际上就是将极地投影的中央经度旋转90度，示例脚本：

1. f = addfile('D:/Temp/nc/seaice_conc_daily_sh_f17_20080830_v02r00.nc')
   
2. data = f['seaice_conc_cdr'][0,:,:]
   
3. #Plot
   
4. mproj = projinfo(proj='stere', lon_0=-90, lat_0=-90, lat_ts=-70)
   
5. axesm(projinfo=mproj, gridline=True)
   
6. lworld = shaperead('D:/Temp/map/country1.shp')
   
7. geoshow(lworld)
   
8. layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)
   
9. colorbar(layer)
   
10. t = f.gettime(0)
    
11. title('Sea ice concentration (' + t.strftime('%Y-%m-%d') + ')')

f = addfile('D:/Temp/nc/seaice_conc_daily_sh_f17_20080830_v02r00.nc')<br /> data = f['seaice_conc_cdr'][0,:,:]<br /> #Plot<br /> mproj = projinfo(proj='stere', lon_0=-90, lat_0=-90, lat_ts=-70)<br /> axesm(projinfo=mproj, gridline=True)<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(lworld)<br /> layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)<br /> colorbar(layer)<br /> t = f.gettime(0)<br /> title('Sea ice concentration (' + t.strftime('%Y-%m-%d') + ')')

在Console中输入f.proj可以显示出从文件中读出的投影信息（Proj4投影字符串，参考此网页：http://remotesensing.org/geotiff/proj_list/）作为参考：
>>> f.proj
+proj=stere +lon_0=0.0 +lat_0=-90.0 +lat_ts=-70.0 +k=1.0 +x_0=0.0 +y_0=0.0

你add point的代码是什么样的？

MeteoInfo

real-slimshady 发表于 2016-4-5 16:08
老师，如何让图像顺时针或逆时针旋转90°？还有我用了add point但是无法添加点

旋转及添加点的代码：

1. f = addfile('D:/Temp/nc/seaice_conc_daily_sh_f17_20080830_v02r00.nc')
   
2. data = f['seaice_conc_cdr'][0,:,:]
   
3. #Plot
   
4. mproj = projinfo(proj='stere', lon_0=-90, lat_0=-90, lat_ts=-70)
   
5. axesm(projinfo=mproj, gridline=True)
   
6. lworld = shaperead('D:/Temp/map/country1.shp')
   
7. geoshow(lworld)
   
8. layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)
   
9. colorbar(layer)
   
10. t = f.gettime(0)
    
11. title('Sea ice concentration (' + t.strftime('%Y-%m-%d') + ')')
    
12. #Add point
    
13. lon = -100
    
14. lat = -80
    
15. geoshow(lat, lon, size=14, color='g', marker='S')

f = addfile('D:/Temp/nc/seaice_conc_daily_sh_f17_20080830_v02r00.nc')<br /> data = f['seaice_conc_cdr'][0,:,:]<br /> #Plot<br /> mproj = projinfo(proj='stere', lon_0=-90, lat_0=-90, lat_ts=-70)<br /> axesm(projinfo=mproj, gridline=True)<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(lworld)<br /> layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)<br /> colorbar(layer)<br /> t = f.gettime(0)<br /> title('Sea ice concentration (' + t.strftime('%Y-%m-%d') + ')')<br /> #Add point<br /> lon = -100<br /> lat = -80<br /> geoshow(lat, lon, size=14, color='g', marker='S')

real-slimshady

MeteoInfo 发表于 2016-4-5 16:52
旋转及添加点的代码：

老师，我用这个程序可以画出日平均数据的图，但对于月平均数据会报错，同时想请问如何求文件中数据的平均值

1. f = addfile('D:/nsidc/seaice_conc_monthly_sh_f17_200801_v02r00.nc')
   
2. data = f['seaice_conc_cdr'][0,:,:]
   
3. #Plot
   
4. mproj = projinfo(proj='stere', lon_0=-90, lat_0=-90, lat_ts=-70)
   
5. axesm(projinfo=f.proj, gridline=True)
   
6. lworld = shaperead('E:/MeteoInfo/map/country1.shp')
   
7. geoshow(lworld)
   
8. layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)
   
9. colorbar(layer)
   
10. t = f.gettime(0)
    
11. title('Sea ice concentration (' + t.strftime('%Y-%m') + ')')
    
12. #Add point
    
13. lon = 76.3
    
14. lat = -69.3
    
15. geoshow(lat, lon, size=14, color='k', marker='S')

复制代码

MeteoInfo

real-slimshady 发表于 2017-5-4 16:24
老师，我用这个程序可以画出日平均数据的图，但对于月平均数据会报错，同时想请问如何求文件中数据的平均 ...

应该是变量名的问题，日均文件中的变量名和月均文件中的变量名可能不同。如果刚接触MeteoInfoLab，建议仔细看看tutorials：http://www.meteothinker.com/docs ... uide/tutorials.html

MeteoInfo

real-slimshady 发表于 2017-5-4 16:24
老师，我用这个程序可以画出日平均数据的图，但对于月平均数据会报错，同时想请问如何求文件中数据的平均 ...

参考此脚本：

1. fn = 'C:/Temp/seaice_conc_monthly_sh_f17_200801_v02r00.nc'
   
2. f = addfile(fn)
   
3. data = f['seaice_conc_monthly_cdr'][0,:,:]
   
4. data[data>1] = nan
   
5. #Plot
   
6. axesm(projinfo=f.proj, gridline=True)
   
7. lworld = shaperead('D:/Temp/map/country1.shp')
   
8. geoshow(lworld, edgecolor='k', facecolor='lightgray')
   
9. layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)
   
10. colorbar(layer)

fn = 'C:/Temp/seaice_conc_monthly_sh_f17_200801_v02r00.nc'<br /> f = addfile(fn)<br /> data = f['seaice_conc_monthly_cdr'][0,:,:]<br /> data[data>1] = nan<br /> #Plot<br /> axesm(projinfo=f.proj, gridline=True)<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(lworld, edgecolor='k', facecolor='lightgray')<br /> layer = imshowm(data, 20, cmap='WhBlGrYeRe', proj=f.proj)<br /> colorbar(layer)

langit

楼主，我按照你的运行，为什么colorbar显示的不是0~1，是到2.5啊大陆上的值也显示出来了 ，是我数据的问题吗？数据下载的是nsidc的每天的海冰密集度。

MeteoInfo

langit 发表于 2017-5-25 11:22
楼主，我按照你的运行，为什么colorbar显示的不是0~1，是到2.5啊大陆上的值也显示出来了 ，是我数 ...

注意楼上的这句代码：
data[data>1] = nan

langit

MeteoInfo 发表于 2017-5-25 11:32
注意楼上的这句代码：
data[data>1] = nan

嗯 我加过这个语句 但是画出来的图 感觉不大对 是这样的 您看一下