MeteoInfoLab脚本示例：IMS冰雪数据

MeteoInfo

示例读取和绘制NSIDC的IMS冰雪数据（ASCII格式，24km分辨率）：

1. #Data source: https://nsidc.org/data/g02156#ascii
   
2. print 'Read data...'
   
3. fn = 'D:/Temp/ascii/ims1997035_24km_v1.1.asc'
   
4. f = open(fn)
   
5. data = zeros(1024 * 1024, dtype='int')
   
6. li = 0
   
7. ii = 0
   
8. for line in f.read().splitlines():
   
9.     if li < 30:    #Head lines
   
10.         print line
    
11.     else:          #Data
    
12.         for d in line:
    
13.             data[ii] = int(d)
    
14.             ii += 1
    
15.     li += 1
    
16. f.close()
    
17. data = data.reshape([1024,1024])
    
18. 
    
19. #Set x, y
    
20. dx = 23684.997
    
21. x = arange1(-12126597.0+dx/2, 1024, dx)
    
22. y = arange1(-12126840.0+dx/2, 1024, dx)
    
23. 
    
24. #Plot
    
25. print 'Plot...'
    
26. proj = projinfo(proj='stere', lat_0=90, lon_0=-80, lat_ts=60)
    
27. axesm(projinfo=proj, frameon=False, griddx=30, griddy=30, gridline=True)
    
28. geoshow('country', edgecolor='k')
    
29. cols = makecolors(5)
    
30. cols[0] = 'w'
    
31. layer = imshowm(x, y, data, colors=cols, proj=proj)
    
32. colorbar(layer, ticklabels=['None','Sea','Land','Snow','Sea ice'])

#Data source: <a href=\"https://nsidc.org/data/g02156#ascii\" target=\"_blank\">https://nsidc.org/data/g02156#ascii</a><br /> print 'Read data...'<br /> fn = 'D:/Temp/ascii/ims1997035_24km_v1.1.asc'<br /> f = open(fn)<br /> data = zeros(1024 * 1024, dtype='int')<br /> li = 0<br /> ii = 0<br /> for line in f.read().splitlines():<br />     if li < 30:    #Head lines<br />         print line<br />     else:          #Data<br />         for d in line:<br />             data[ii] = int(d)<br />             ii += 1<br />     li += 1<br /> f.close()<br /> data = data.reshape([1024,1024])<br /> <br /> #Set x, y<br /> dx = 23684.997<br /> x = arange1(-12126597.0+dx/2, 1024, dx)<br /> y = arange1(-12126840.0+dx/2, 1024, dx)<br /> <br /> #Plot<br /> print 'Plot...'<br /> proj = projinfo(proj='stere', lat_0=90, lon_0=-80, lat_ts=60)<br /> axesm(projinfo=proj, frameon=False, griddx=30, griddy=30, gridline=True)<br /> geoshow('country', edgecolor='k')<br /> cols = makecolors(5)<br /> cols[0] = 'w'<br /> layer = imshowm(x, y, data, colors=cols, proj=proj)<br /> colorbar(layer, ticklabels=['None','Sea','Land','Snow','Sea ice'])

chongzika

120w进到圈里了，不大美观

hello-lyn

感谢分享，最近在处理这个数据，另外想请教，如何提取出某一个区域（已知经纬度范围）的IMS雪盖值呢？