标签：
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)
标签：
原文地址：http://www.cnblogs.com/yaqiang/p/4601949.html
