MeteoInfoLab脚本示例：inpolygon

MeteoInfo

inpollygon函数是用来判断带坐标（x/y）的数据是否在某个或者一组多边形（Polygon）中，返回的结果中如果做多边形内则值为1，否则值为-1。下面一个例子演示了利用一个shape文件和inpolygon函数生成这种0、1数据。

1、生成格点数据的坐标（x/y），arange1函数生成一个1维数组，三个参数分别为起始值，个数和步长。用meshgrid函数生成x/y二维数组。

2、读取shape文件，返回一个图层m_china。

3、利用数组a的inpolygon函数生成一个新数组b，b中多边形内的点的值为True，其它为False。

4、将b转为 int 类型（True -> 1, False -> 0）。

5、利用数组b的savegrid方法将数组b及其坐标信息保存到一个Surfer ASCII grid格式的文件中。

6、绘图，这是为了测试一下。数组b中的1用红色的点表示，0用灰色的点显示。

脚本程序：

1. #Create x/y vectorx = arange1(60, 80, 1)
   
2. y = arange1(10, 50, 1)
   
3. mx, my = meshgrid(x, y)
   
4. 
   
5. #Read shape file
   
6. m_china = shaperead('D:/Temp/map/china.shp')
   
7. 
   
8. #inpolygon function
   
9. b = inpolygon(mx, my, m_china)
   
10. 
    
11. #Convert boolean array to int
    
12. b = b.astype('int')
    
13. 
    
14. #Save grid data
    
15. #b.savegrid(x, y, 'D:/Temp/test.dat')
    
16. 
    
17. #Plot
    
18. axesm()
    
19. geoshow(m_china)
    
20. ss = makesymbolspec('point', {'value':0, 'color':'lightgray', 'size':2, 'edge':False}, {'value':1, 'color':'red', 'size':2, 'edge':False})
    
21. layer = scatterm(x, y, b, symbolspec=ss)

#Create x/y vectorx = arange1(60, 80, 1)<br /> y = arange1(10, 50, 1)<br /> mx, my = meshgrid(x, y)<br /> <br /> #Read shape file<br /> m_china = shaperead('D:/Temp/map/china.shp')<br /> <br /> #inpolygon function<br /> b = inpolygon(mx, my, m_china)<br /> <br /> #Convert boolean array to int<br /> b = b.astype('int')<br /> <br /> #Save grid data<br /> #b.savegrid(x, y, 'D:/Temp/test.dat')<br /> <br /> #Plot<br /> axesm()<br /> geoshow(m_china)<br /> ss = makesymbolspec('point', {'value':0, 'color':'lightgray', 'size':2, 'edge':False}, {'value':1, 'color':'red', 'size':2, 'edge':False})<br /> layer = scatterm(x, y, b, symbolspec=ss)

kmdqlb

王老师好 ，请问这些函数在MeteoInfo中能用吗？谢谢

MeteoInfo

kmdqlb 发表于 2015-6-27 14:37
王老师好 ，请问这些函数在MeteoInfo中能用吗？谢谢

这些是专门为MeteoInfoLab设计的函数，当然MeteoInfo里也都有相应的功能，但实现起来比较麻烦。这也是为什么要开发MeteoInfoLab的一个重要原因。

liusinuo

MeteoInfo 发表于 2015-6-27 16:55
这些是专门为MeteoInfoLab设计的函数，当然MeteoInfo里也都有相应的功能，但实现起来比较麻烦。这也是为 ...

王老师您好，您说的“meteoinfo里也有相应的功能，但实现起来比较麻烦”，您能给我说下大致的思路么？我对meteoinfolab不太熟啊

liusinuo

MeteoInfo 发表于 2015-6-27 16:55
这些是专门为MeteoInfoLab设计的函数，当然MeteoInfo里也都有相应的功能，但实现起来比较麻烦。这也是为 ...

王老师您好，我根据您的代码，只是更改了文件路径，但是出现了下面的错误，请问这是什么原因啊？

qxsw2016

利用此方法，计算某一区域内，最大值，最小值，平均值，脚本仅供参考：
f = addfile('F:/DataSet/ncep/hgt.mon.mean.nc')
var = f['hgt'][0,'500','10:60','60:140']
temp = var
#print var.dimvalue(1)
#x = arange(60, 140, 1.0)
#y = arange(10, 60, 1.0)
x = var.dimvalue(1)
y = var.dimvalue(0)
mx, my = meshgrid(x,y)
#Read shape file
m_china = shaperead('F:/MeteoInfo/MeteoInfo/Map/china.shp')
#inpolygon function
b = inpolygon(mx, my, m_china)
#Convert boolean array to int
b = b.astype('int')
#print b.shape
#Save grid data
#b.savegrid(x, y, 'D:/Temp/test.dat')

#Plot
axesm()
geoshow(m_china)
ss = makesymbolspec('point', {'value':0, 'color':'lightgray', 'size':2, 'edge':False}, {'value':1, 'color':'red', 'size':2, 'edge':False})
#ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6}, \
#    {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6})
layer = scatterm(mx, my,b,symbolspec=ss)
layer = contourm(var)
#layer = contourfm(x, y, b, 20)
colorbar(layer)
nx = var.dimlen(1)
ny = var.dimlen(0)
for j in arange(ny):
    for i in arange(nx):
        if(b[j,i]==0):
            temp[j,i] = nan
#print temp[16,:]
xid = argmax(temp)%33
yid = argmax(temp)/33
#print xid,yid
print 'maxvalue: ',max(temp),'[x,y]',x[xid],y[yid]
xid = argmin(temp)%33
yid = argmin(temp)/33
print 'minvalue: ',min(temp),'[x,y]',x[xid],y[yid]
print 'meanvalue: ',mean(temp)