Записки программиста, обо всем и ни о чем. Но, наверное, больше профессионального.

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);
    }
}


А теперь таки про 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

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

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

Архив блога

Ярлыки

linux (241) python (191) citation (186) web-develop (170) gov.ru (159) video (124) бытовуха (115) sysadm (100) GIS (97) Zope(Plone) (88) бурчалки (84) Book (83) programming (82) грабли (77) Fun (76) development (73) windsurfing (72) Microsoft (64) hiload (62) internet provider (57) opensource (57) security (57) опыт (55) movie (52) Wisdom (51) ML (47) driving (45) hardware (45) language (45) money (42) JS (41) curse (40) bigdata (39) DBMS (38) ArcGIS (34) history (31) PDA (30) howto (30) holyday (29) Google (27) Oracle (27) tourism (27) virtbox (27) health (26) vacation (24) AI (23) Autodesk (23) SQL (23) Java (22) humor (22) knowledge (22) translate (20) CSS (19) cheatsheet (19) hack (19) Apache (16) Manager (15) web-browser (15) Никонов (15) functional programming (14) happiness (14) music (14) todo (14) Klaipeda (13) PHP (13) course (13) scala (13) weapon (13) HTTP. Apache (12) SSH (12) frameworks (12) hero (12) im (12) settings (12) HTML (11) SciTE (11) USA (11) crypto (11) game (11) map (11) HTTPD (9) ODF (9) Photo (9) купи/продай (9) benchmark (8) documentation (8) 3D (7) CS (7) DNS (7) NoSQL (7) cloud (7) django (7) gun (7) matroska (7) telephony (7) Microsoft Office (6) VCS (6) bluetooth (6) pidgin (6) proxy (6) Donald Knuth (5) ETL (5) NVIDIA (5) Palanga (5) REST (5) bash (5) flash (5) keyboard (5) price (5) samba (5) CGI (4) LISP (4) RoR (4) cache (4) car (4) display (4) holywar (4) nginx (4) pistol (4) spark (4) xml (4) Лебедев (4) IDE (3) IE8 (3) J2EE (3) NTFS (3) RDP (3) holiday (3) mount (3) Гоблин (3) кухня (3) урюк (3) AMQP (2) ERP (2) IE7 (2) NAS (2) Naudoc (2) PDF (2) address (2) air (2) british (2) coffee (2) fitness (2) font (2) ftp (2) fuckup (2) messaging (2) notify (2) sharepoint (2) ssl/tls (2) stardict (2) tests (2) tunnel (2) udev (2) APT (1) CRUD (1) Canyonlands (1) Cyprus (1) DVDShrink (1) Jabber (1) K9Copy (1) Matlab (1) Portugal (1) VBA (1) WD My Book (1) autoit (1) bike (1) cannabis (1) chat (1) concurrent (1) dbf (1) ext4 (1) idioten (1) join (1) krusader (1) license (1) life (1) migration (1) mindmap (1) navitel (1) pneumatic weapon (1) quiz (1) regexp (1) robot (1) science (1) serialization (1) spatial (1) tie (1) vim (1) Науру (1) крысы (1) налоги (1) пианино (1)