VSnake notes: spatialReference

2012-10-10

spatialReference

- Где я нахожусь?
- В корзине воздушного шара.

Как и обещал вчера, поговорим о том, куда девается определение системы координат (spatialReference) при передаче FeatureSet в скрипт геопроцессора ArcGIS.

Но сначала закроем вчерашнюю тему, про дырку от бублика. Я нашел только один способ автоматически превращать всякую фигню, нарисованную пользователем, в корректные, с точки зрения ArcGIS, фигуры. Этот способ заключается в использовании Geometry Service «Simplify». И не спрашивайте, почему он так называется - «упростить». Он делает ровно то, что мне нужно — превращает дырки в нормальные полигоны и делает из самопересекающегося полигона правильный составной полигон.

Чтобы использовать этот сервис в прикладухе на ArcGIS Silverlight API, достаточно написать нечто вроде

private void evtDrawComplete(object sender, DrawEventArgs e) {
    currGeom = e.Geometry;
    draw.DrawMode = DrawMode.None;
    draw.IsEnabled = false;
    try {
        var gl = getRLLayer();
        var graphic = new ESRI.ArcGIS.Client.Graphic() {
            Geometry = currGeom, Symbol = polySymbol
        };
        gl.Graphics.Add(graphic);

        var geometryService = new GeometryService("http://tasks.arcgisonline.com/ArcGIS/rest/services/"+
            "Geometry/GeometryServer");
        var graphicList = new List<Graphic>();
        graphicList.Add(graphic);

        geometryService.SimplifyCompleted += (sndr, args) => {
            gl.Graphics.Remove(graphic);
            graphic.Geometry = args.Results[0].Geometry;
            gl.Graphics.Add(graphic);
            askGeoprocessor(graphic.Geometry); // send good polygon to geoprocessor
        };

        geometryService.Failed += (sndr, args) => {
            MessageBox.Show("Сбой нормализации полигона: " + args.Error);
        };

        geometryService.SimplifyAsync(graphicList);
    }
    catch(Exception ex) {
        string msg = string.Format("Сбой отправки запроса: \n [{0}]", ex.Message);
        MessageBox.Show(msg);
    }
}

http://help.arcgis.com/en/webapi/silverlight/samples/start.htm#Simplify

http://tasks.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer/simplify

А теперь таки про spatialReference.

Фишка в том, что в Python скрипте геопроцессора нет доступа к сведениям о системе координат переданного в скрипт полигона (используется FeatureSet), хотя сведения эти передаются явным образом

inputPolygon={
"geometryType":"esriGeometryPolygon",
"spatialReference":{"wkid":102100},
"features":[{
 "geometry":{
  "spatialReference":{"wkid":102100},
  "rings":[[[7592337.47835702,9803507.48815798],[7924991.42545401,10312272.348424],[8277213.25179201,9979618.40132698],[7592337.47835702,9803507.48815798]]]
 },
 "attributes":{}
}]}

Проведенные тесты показывают, что где-то глубоко внутре используемых обьектов эта информация есть, равно как и координаты дырки от бублика. Просто через API не выходит получить прямой доступ к данным.

Но важно не это, важно другое — что же делать, трахтибидох? Какой смысл в координатах если мы не знаем систему координат, точку отсчета?

На мой взгляд, сделать можно две вещи: передавать spatialreference отдельным параметром и/или трансформировать переданные координаты в известную и желаемую систему координат. Второй способ лично мне нравится больше.

Вот смотрите

import arcpy
from arcpy import env

fsetObj = arcpy.GetParameter(0)
fsetDesc = arcpy.Describe(fsetObj)
if hasattr(fsetDesc, 'spatialReference'): log.info("arcpyStuff, input fset spatialReference '%s'" % fsetDesc.spatialReference) # haven't

WKID = 4326 # WGS-1984 http://anothergisblog.blogspot.com/2011/05/spatial-reference-class.html
sr = arcpy.SpatialReference()
sr.factoryCode = WKID
sr.create()
# env.outputCoordinateSystem = sr

# или возьмем spatialreference из заранее подготовленного featureclass
gdbPath = r'''\\cache\MXD\seismo\Seis_button.gdb'''
seisFCName = r'''APP_GP_SEISM2D_L'''
seisDesc = arcpy.Describe(os.path.join(gdbPath, seisFCName))
log.info("arcpyStuff, seismoprofiles WKID '%s'" % (seisDesc.spatialReference.factoryCode))

rows = arcpy.SearchCursor(fsetObj, '', sr) # decimal degree
# rows = arcpy.SearchCursor(fsetObj, '', seisDesc.spatialReference) # meters
for row in rows:
    geom = row.shape
    log.info("arcpyStuff searchcursor, geometry geojson '%s'" % (geom.__geo_interface__)) # {'type': 'Polygon', 'coordinates': []}
    if hasattr(geom, 'spatialReference'): log.info("arcpyStuff searchcursor, geometry spatialReference '%s'" % (geom.spatialReference)) # haven't
    if geom.partCount <= 0 or geom.area <=0:
        raise NameError("Wrong input polygon, you should send no selfintersected clockwise drawed single ring")
    for part in geom:
        log.info("arcpyStuff searchcursor, geometry part, len '%s'" % (len(part)))
        for pnt in part:
            log.info("arcpyStuff searchcursor, point x '%s', y '%s'" % (pnt.X, pnt.Y))

Трансформация происходит (задается) на строке

rows = arcpy.SearchCursor(fsetObj, '', sr) # decimal degree

И, как ни странно, это работает. На входе имеем

WKID = 102100
geometry geojson '{'type': 'Polygon', 'coordinates': [[(7924991.4255000018, 10390543.865400001), (8257645.3726000004, 9803507.4882000014), (7631473.2368000001, 9901346.8844000027), (7924991.4255000018, 10390543.865400001)]]}'

а на выходе

WKID = 4326
geometry geojson '{'type': 'Polygon', 'coordinates': [[(71.191409240428044, 67.809244466732096), (74.17969049045422, 65.730626450487776), (68.554690489665163, 66.089364228422752), (71.191409240428044, 67.809244466732096)]]}'

Легко проверить, что это одни и те же координаты, только в разных СК.

original post http://vasnake.blogspot.com/2012/10/spatialreference.html

Tools

VSnake notes

2012-10-10

spatialReference

Комментариев нет:

Отправить комментарий

Архив блога

Ярлыки

Обо мне

Links