MeteoInfoLab脚本示例：inpolygon

inpollygon函数是用来判断带坐标（x/y）的数据是否在某个或者一组多边形（Polygon）中，返回的结果中如果做多边形内则值为1，否则值为-1。下面一个例子演示了利用一个shape文件和inpolygon函数生成这种0、1数据。需要下载最新的MeteoInfo版本（1.2.6R1）：http://yun.baidu.com/share/link?shareid=669776748&uk=51062435

1、生成一个格点数据（2维数组），可以用zeros函数：
a
= zeros((50,80),
dtype=‘int‘)
生成的数组a是50*80的2维数组，其值均为0。

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

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

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

5、将b中-1改变为0。

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

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

脚本程序：

#Create a 50*80 2D array
a = zeros((50,80), dtype=‘int‘)
#Create x/y vector
x = arange1(60, 80, 1)
y = arange1(10, 50, 1)
#Read shape file
m_china = shaperead(‘D:/Temp/map/china.shp‘)
#inpolygon function, in: 1, out: -1
b = a.inpolygon(x, y, m_china)
#Change -1 to 0
b[b==-1] = 0
#Save to a surfer ASCII grid data file
fn = ‘D:/Temp/test/sdata.dat‘
b.savegrid(x, y, fn)
#Plot
axesm()
geoshow(m_china)
ss = makesymbolspec(‘point‘, {‘value‘:0, ‘color‘:‘lightgray‘, ‘size‘:2, ‘edge‘:False}, {‘value‘:1, ‘color‘:‘red‘, ‘size‘:2, ‘edge‘:False})
layer = scatterm(x, y, b, symbolspec=ss)