Le portail francophone de la géomatique

from osgeo import osr

#Dictionnaires python décrivant des projections de référence

Lambert2E={
"Name":"Lambert 2 étendu",
"Datum":"NTF",
"Spheroid":"Clarke1880",#+ellps=clrk80
"Sph_SemiMajor":6378249.2,#+a
"Sph_SemiMinor":6356515,#+b
"Sph_Inverse_Flat":293.466,#+rf
"Sph_Prime_Meridian":2.337229167,#Paris, +pm
"Sph_Angular_Unit_Name":"Degree",
"Sph_Angular_Unit_Scale":0.0174532925199433,
"Projection":"Lambert_Conformal_Conic",#+proj=lcc
"Prj_False_Easting":600000,#+x_0
"Prj_False_Northing":2200000,#+y_0
"Prj_Central_Meridian":0,#+long_0
"Prj_Standard_Parallel_1":45.8989188889,#+lat_1
"Prj_Standard_Parallel_2":47.6960144444,#+lat_2
"Prj_Latitude_Of_Origin":46.8,#+lat_0
"Prj_Scale_Factor":1,#+k_0
"Prj_Unit_Name":"Meter",#+units
"Prj_Unit_Scale":1#+units
}

Lambert93={
"Name":"RGF93_Lambert_93",
"Datum":"RGF93",
"Spheroid":"GRS_1980",
"Sph_SemiMajor":6378137.0,
"Sph_SemiMinor":6356752.31414,
"Sph_Inverse_Flat":298.257222101,
"Sph_Prime_Meridian":0,#Greenwitch
"Sph_Angular_Unit_Name":"Degree",
"Sph_Angular_Unit_Scale":0.0174532925199433,
"Projection":"Lambert_Conformal_Conic",
"Prj_False_Easting":700000,
"Prj_False_Northing":6600000,
"Prj_Central_Meridian":3,
"Prj_Standard_Parallel_1":44,
"Prj_Standard_Parallel_2":49,
"Prj_Latitude_Of_Origin":46.5,
"Prj_Scale_Factor":1,
"Prj_Unit_Name":"Meter",
"Prj_Unit_Scale":1
}

#fonction la + intéressante, parse à l'aide d'osr le fichier prj et retourne un dictionnaire formaté comme ci-dessus
def parsePrj(projWKT):
    srs = osr.SpatialReference(projWKT)
    proj={}
    proj["Name"]=srs.GetAttrValue("PROJCS")
    proj["Datum"]=srs.GetAttrValue("DATUM")
    proj["Spheroid"]=srs.GetAttrValue("SPHEROID")
    proj["Sph_SemiMajor"]=srs.GetSemiMajor()#srs.GetAttrValue("SPHEROID", 1)
    proj["Sph_SemiMinor"]=srs.GetSemiMinor()
    proj["Sph_Inverse_Flat"]=float(srs.GetAttrValue("SPHEROID", 2))
    proj["Sph_Prime_Meridian"]=float(srs.GetAttrValue("PRIMEM",1))
    proj["Sph_Angular_Unit_Name"]=srs.GetAttrValue("GEOGCS|UNIT")
    proj["Sph_Angular_Unit_Scale"]=srs.GetAngularUnits()
    proj["Projection"]=srs.GetAttrValue("PROJECTION")
    proj["Prj_False_Easting"]=srs.GetProjParm("false_easting")
    proj["Prj_False_Northing"]=srs.GetProjParm("false_northing")
    proj["Prj_Central_Meridian"]=srs.GetProjParm("central_meridian")
    proj["Prj_Standard_Parallel_1"]=srs.GetProjParm("standard_parallel_1")
    proj["Prj_Standard_Parallel_2"]=srs.GetProjParm("standard_parallel_2")
    proj["Prj_Latitude_Of_Origin"]=srs.GetProjParm("latitude_of_origin")
    proj["Prj_Scale_Factor"]=srs.GetProjParm("Scale_Factor")
    proj["Prj_Unit_Name"]=srs.GetLinearUnitsName()
    proj["Prj_Unit_Scale"]=srs.GetLinearUnits()
    return proj

#comparer quelques paramètres pr déterminer si proj1 == proj2
def compareProjSys(proj1, proj2):
    return proj1["Prj_False_Easting"] == proj2["Prj_False_Easting"] and proj1["Prj_False_Northing"] == proj2["Prj_False_Northing"]

#fonction principale pour déterminer la projection
#l'objectif ici était surtout de déterminer si c'est du lambert 2e
#les autres projections n'étant pas traité par la suite script, il n'était pas nécessaire de les déterminer précisément
def estimateProjSys(prjFile):
    lb2eNameVariantes=["NTF_Lambert_II_étendu", "NTF_Lambert_II_Carto", "Lambert II Carto", "NTF_France_II_degrees"]#trouvé empiriquement
    try:
        prjFile=open(prj,"r")
        projWKT=prjFile.read()
        prjFile.close()
    except:
        return "prj err"
    for variante in lb2eNameVariantes:
        if variante in projWKT:
            return "NTF Lb2e"
    #on check si c'est pas autre chose que du NTF
    if "RGF" in projWKT:
        return "RGF93"#n'importe quel system basé sur RGF93...
    elif "WGS" in projWKT:
        return "WGS84"
    else:#Pas de nom de projection ou nom de datum pour nous aiguiller, on parse le fichier pour ensuite comparer certain paramètres
        unknowProj=parsePrj(projWKT)
        if compareProjSys(Lambert2E, unknowProj):
            return "NTF Lb2e (checked)"
        elif "NTF" in projWKT:
            return "NTF unknow"
        elif compareProjSys(Lambert93, unknowProj):
            return "RGF93 (checked)"
        else:
            return "Other"