- 积分
- 55950
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-6-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
借用坛友的问题:如何判断点是否在某区域内?这里提供一个更完整的脚本程序来生成一个等间隔网格点图层,判断每个点是否在青藏高原内并在点图层中赋相应的属性('Y' or 'N'),输出内部的点的坐标,设置点图层的LegendScheme(区域内的点红色,外面的点绿色)。可以很直观的验证判断结果。
运行结果:
脚本代码:
- #---- Import clr and classes
- import clr
- clr.AddReferenceByPartialName("System")
- clr.AddReferenceByPartialName("System.Drawing")
- from System import *
- from System.Drawing import *
- clr.AddReference("MeteoInfoC.dll")
- from MeteoInfoC import *
- from MeteoInfoC.Data.MapData import *
- from MeteoInfoC.Geoprocess import *
- from MeteoInfoC.Layer import *
- from MeteoInfoC.Shape import *
- from MeteoInfoC.Legend import *
- #---- Open layer
- mapdir = "D:\\MyProgram\\CSharp\\MeteoInfoDev\\MeteoInfo\\bin\\Debug\\Map\"
- aFile = mapdir + "tibet_ASCII.wmp"
- aLayer = MapDataManage.OpenLayer(aFile)
- #---- Add layer
- mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer)
- mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer.Handle, 0)
- #---- Create a grid point and determine grid points are in or out of Tibet
- nLayer = VectorLayer(ShapeTypes.Point)
- nLayer.LayerName = "New_Point_Layer"
- nLayer.Visible = True
- fieldName = "Inside"
- nLayer.EditAddField(fieldName, String)
- print 'The points inside Tibet are:'
- tibetArea = aLayer.ShapeList[0]
- sLon = 62.0
- sLat = 25.0
- delt = 1.0
- for i in range(0, 45):
- lon = sLon + i * delt
- for j in range(0, 24):
- aPS = PointShape()
- lat = sLat + j * delt
- aPoint = PointD(lon, lat)
- aPS.Point = aPoint
- shapeNum = nLayer.ShapeNum
- print "Shape Number: %d" % (shapeNum)
- nLayer.EditInsertShape(aPS, shapeNum)
- if GeoComputation.PointInPolygon(tibetArea, aPoint):
- nLayer.EditCellValue(fieldName, shapeNum, "Y")
- print "Lon: %5.1f, Lat: %5.1f" % (lon, lat)
- else:
- nLayer.EditCellValue(fieldName, shapeNum, "N")
- aLS = nLayer.CreateLegendScheme(LegendType.UniqueValue, fieldName)
- aLS.LegendBreaks[0].DrawOutline = False
- aLS.LegendBreaks[0].Color = Color.Green
- aLS.LegendBreaks[1].DrawOutline = False
- aLS.LegendBreaks[1].Color = Color.Red
- nLayer.LegendScheme = aLS
- print "Layer Name: %s" % (nLayer.LayerName)
- mipy.MapDocument.ActiveMapFrame.AddLayer(nLayer)
- print 'Finished!'
复制代码
|
|