2017-05-02 2 views
0

J'ai une liste de géocroords (lat, long) et un shapefile avec différentes couches. Je veux être capable d'identifier à quelle couche appartient chaque coord.Conversion de la latitude du shapefile en vrai

Mais le shapefile (SHP) a poligons dans lequel la latitude et la longitude sont exprimées en nombre dans une plage bizarre, comme 120724,86008864 et 484497,34058312

Je sais que le fichier PRJ contient les informations sur la façon dont cette transformation a été fait, mais je ne semble pas obtenir comment. Voici:

PROJCS [ "RD_New", GEOGCS [ "GCS_Amersfoort", DONNÉE [ "D_Amersfoort", SPHEROID [ "Bessel_1841", 6377397.155,299.1528128]], PRIMEM [ "Greenwich", 0], UNITÉ [ » Degré ", 0.0174532925199432955]], PROJECTION [" Double_Stereographic "], PARAMETRE [" False_East ", 155000], PARAMETRE [" False_Northing ", 463000], PARAMETRE [" Central_Meridian ", 5.38763888888889], PARAMETRE [" Scale_Factor ", 0.9999079], PARAMETRE ["Latitude_Of_Origin", 52.15616055555555], UNIT ["Compteur", 1]]

La question concrète est de savoir comment transformer un point lat/long régulier en ceux du fichier de formes.

Travailler en Python, avec cette bibliothèque http://gdal.org/python/

Toute aide est appréciée.

+0

Vous voulez obtenir les coordonnées? Ou analyser le fichier prj? –

+0

Je peux analyser le fichier sans problème, mais les lat/long dans le fichier de formes sont dans une échelle ou une plage différente. Par exemple, une latitude de 120724.86008864 qui est manifestement erronée. J'ai des latitudes normales telles que 52.3605883. Donc je veux savoir quelle transformation je dois appliquer à mes coordonnées normales pour être comme celles du fichier. Ensuite, je serai en mesure de les croiser avec les poligones dans les couches. –

Répondre

0
# define input 
shape_file = "file.shp" 
o_lat = 52.3605883 
o_lon = 4.8593157 

# geospatial bureocracy 
driver = ogr.GetDriverByName('ESRI Shapefile') 
shape = driver.Open(shape_file) 
layer = shape.GetLayer() 
geo_ref = layer.GetSpatialRef() 
point_ref = ogr.osr.SpatialReference() 
point_ref.ImportFromEPSG(4326) 
ctran = ogr.osr.CoordinateTransformation(point_ref, geo_ref) 

# critical part: transform longitude/latitude to the shapefile's projection 
[t_lon, t_lat, z] = ctran.TransformPoint(o_lon, o_lat) 
print('original coords', o_lon, o_lat) 
print('transformed coords', t_lon, t_lat) 

# create the needle 
point = ogr.Geometry(ogr.wkbPoint) 
point.SetPoint_2D(0, t_lon, t_lat) 
layer.SetSpatialFilter(point) 

# look it up 
for feature in layer: 
    polygon = feature.GetGeometryRef() 
    if polygon.Contains(point): 
     print('Found it', feature.ExportToJson() 
0
 import re 
     regex = "\d{1,3}\.\d+" 
     s ="""PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199432955]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",155000],PARAMETER["False_Northing",463000],PARAMETER["Central_Meridian",5.38763888888889],PARAMETER["Scale_Factor",0.9999079],PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1]] """ 

     m = re.search(regex, s) 

     if m: 
      print m.groups() 
+0

Quelle est votre intention? –

+0

Je veux obtenir tout décimal avec max 3 chiffres après le point décimal de la chaîne. ex: 12,55 3,44 321,11 –