Improved bugs and added proj/datu/ellip names leading to much better test coverage (70-90% of the crs test cases)

pycrs/__init__.py

@@ -1,3 +1,11 @@
+"""
+PyCRS
+Karim Bahgat, 2015
+
+Intro...
+"""
+
+
 from . import loader
 from . import parser
 from . import utils

pycrs/elements/containers.py

@@ -1,6 +1,7 @@
 
 from . import directions
 from . import datums
+from . import ellipsoids
     
 # the final CRS object which is instantiated with all of the below and parameters
 # remember to use +no_defs when outputting to proj4
@@ -61,6 +62,9 @@ class Datum:
             return "%s %s" % (self.ellips.to_proj4(), self.datumshift.to_proj4())
         elif isinstance(self.name, datums.Unknown):
             return "%s" % self.ellips.to_proj4()
+        elif not self.name.proj4:
+            # has no proj4 equivaent and is better left unspecified, so only return ellips
+            return "%s" % self.ellips.to_proj4()
         else:
             return "+datum=%s %s" % (self.name.proj4, self.ellips.to_proj4())
 
@@ -104,7 +108,14 @@ class Ellipsoid:
         self.inv_flat = inv_flat
 
     def to_proj4(self):
-        return "+ellps=%s +a=%s +f=%s" % (self.name.proj4, self.semimaj_ax, self.inv_flat)
+        if isinstance(self.name, ellipsoids.Unknown):
+            # has no proj4 equivaent and is better left unspecified
+            return "+a=%s +f=%s" % (self.semimaj_ax, self.inv_flat)
+        elif not self.name.proj4:
+            # has no proj4 equivaent and is better left unspecified
+            return "+a=%s +f=%s" % (self.semimaj_ax, self.inv_flat)
+        else:
+            return "+ellps=%s +a=%s +f=%s" % (self.name.proj4, self.semimaj_ax, self.inv_flat)
 
     def to_ogc_wkt(self):
         return 'SPHEROID["%s", %s, %s]' % (self.name.ogc_wkt, self.semimaj_ax, self.inv_flat)

pycrs/elements/datums.py

@@ -58,13 +58,13 @@ class SphereArcInfo:
     ogc_wkt = "D_Sphere_ARC_INFO" # confirmed but odd that uses D_
     esri_wkt = "D_Sphere_ARC_INFO"
 
-    semimaj_ax = 6370997.0
-    inv_flat = 0.0
+    ellipsdef = "" # ellipsoids...
+    to_wgs84 = None
 
 class Unknown:
-    proj4 = "unknown" # no datum name, just ellips + towgs84 params...
+    proj4 = "" # no datum name, just ellips + towgs84 params...
     ogc_wkt = "Unknown"
     esri_wkt = "Unknown"
 
     ellipsdef = "" # ellipsoids.WGS72()
-    to_wgs84 = None,None
+    to_wgs84 = None

pycrs/elements/ellipsoids.py

@@ -79,7 +79,30 @@ class SphereArcInfo:
     semimaj_ax = 6370997.0
     inv_flat = 0.0
 
+class Krassowsky1940:
+    proj4 = "krass"
+    ogc_wkt = "Krassowsky 1940"
+    esri_wkt = "Krassowsky_1940"
 
+    semimaj_ax = 6378245.0
+    inv_flat = 298.3
+
+class Bessel1841:
+    proj4 = "bessel"
+    ogc_wkt = "Bessel 1841"
+    esri_wkt = "Bessel_1841"
+
+    semimaj_ax = 6377397.155
+    inv_flat = 299.1528128
+
+class Unknown:
+    proj4 = ""
+    ogc_wkt = "Unknown"
+    esri_wkt = "Unknown"
+
+    # values have to be set manually in Ellipsoid class
+    semimaj_ax = None
+    inv_flat = None
 
 
 

pycrs/elements/parameters.py

@@ -24,6 +24,7 @@
 
 from . import datums
 from . import directions
+from . import units
 
 ################
 
@@ -333,7 +334,10 @@ class Unit:
         self.metermultiplier = metermultiplier
 
     def to_proj4(self):
-        return "%s %s" %(self.unittype.to_proj4(), self.metermultiplier.to_proj4())
+        if isinstance(self.unittype, units.Unknown):
+            return "%s" % self.metermultiplier.to_proj4()
+        else:
+            return "%s %s" %(self.unittype.to_proj4(), self.metermultiplier.to_proj4())
 
     def to_ogc_wkt(self):
         return 'UNIT["%s", %s]' %(self.unittype.to_ogc_wkt(), self.metermultiplier.to_ogc_wkt())

pycrs/elements/projections.py

@@ -34,8 +34,8 @@ class UTM:
 
 class ObliqueMercator:
     proj4 = "omerc"
-    ogc_wkt = "Hotine_Oblique_Mercator"
-    esri_wkt = "Hotine_Oblique_Mercator_Azimuth_Natural_Origin"
+    ogc_wkt = "Hotine_Oblique_Mercator_Two_Point_Natural_Origin" #"Hotine_Oblique_Mercator"
+    esri_wkt = "Hotine_Oblique_Mercator_Two_Point_Natural_Origin" #"Hotine_Oblique_Mercator_Azimuth_Natural_Origin"
     
 class AlbersEqualArea:
     proj4 = "aea"
@@ -49,8 +49,20 @@ class CylindricalEqualArea:
     
 class EquiDistantConic:
     proj4 = "eqdc"
-    ogc_wkt = "Equidistant_Conic"
-    esri_wkt = "Equidistant_Conic"
+    ogc_wkt = "Equidistant_Conic" 
+    esri_wkt = "Equidistant_Conic" 
+
+class EquiDistantCylindrical:
+    # same as equirectangular...?
+    proj4 = "eqc"
+    ogc_wkt = "Equidistant_Cylindrical"
+    esri_wkt = "Equidistant_Cylindrical"
+
+class EquiRectangular:
+    # same as equidistant cylindrical
+    proj4 = "eqc"
+    ogc_wkt = "Equirectangular"
+    esri_wkt = "Equirectangular"
 
 class TransverseMercator:
     proj4 = "tmerc"
@@ -85,7 +97,7 @@ class Mollweide:
 class ObliqueStereographic:
     proj4 = "sterea"
     ogc_wkt = "Oblique_Stereographic"
-    esri_wkt = "Stereographic_North_Pole"
+    esri_wkt = "Oblique Stereographic" #"Stereographic_North_Pole"
 
 class Orthographic:
     proj4 = "ortho"
@@ -112,11 +124,6 @@ class VanDerGrinten:
     ogc_wkt = "VanDerGrinten"
     esri_wkt = "Van_der_Grinten_I"
 
-class Equirectangular:
-    proj4 = "eqc"
-    ogc_wkt = "Equirectangular"
-    esri_wkt = "Equidistant_Cylindrical"
-
 class LambertConformalConic:
     proj4 = "lcc"
     ogc_wkt = "Lambert_Conformal_Conic" # possible has some variants
@@ -162,8 +169,30 @@ class ObliqueCylindricalEqualArea:
     ogc_wkt = "Oblique_Cylindrical_Equal_Area"
     esri_wkt = "Oblique_Cylindrical_Equal_Area"
 
-        
+class Polyconic:
+    proj4 = "poly"
+    ogc_wkt = "Polyconic"
+    esri_wkt = "Polyconic"
 
+class EckertIV:
+    proj4 = "eck4"
+    ogc_wkt = "Eckert_IV"
+    esri_wkt = "Eckert_IV"
+
+class EckertVI:
+    proj4 = "eck6"
+    ogc_wkt = "Eckert_VI"
+    esri_wkt = "Eckert_VI"
+
+class AzimuthalEquidistant:
+    proj4 = "aeqd"
+    ogc_wkt = "Azimuthal_Equidistant"
+    esri_wkt = "Azimuthal_Equidistant"
+
+class GeostationarySatellite:
+    proj4 = "geos"
+    ogc_wkt = "Geostationary_Satellite"
+    esri_wkt = "Geostationary_Satellite"
 
 
 

pycrs/elements/units.py

@@ -15,7 +15,7 @@ def find(unitname, crstype, strict=False):
                 if unitname == itemname:
                     return item
                 # special handling of wkt meters which has multiple possibilities
-                elif isinstance(item(), Meter) and crstype.endswith("wkt") and not strict and unitname in ("meters","meter","metre"):
+                elif isinstance(item(), Meter) and crstype.endswith("wkt") and not strict and unitname in ("meters","meter","metre","m"):
                     return item
         except:
             pass
@@ -39,3 +39,8 @@ class Feet:
     proj4 = "..."
     ogc_wkt = "Foot_US"
     esri_wkt = "Foot_US"
+
+class Unknown:
+    proj4 = ""
+    ogc_wkt = "Unknown"
+    esri_wkt = "Unknown"

pycrs/loader.py

@@ -1,3 +1,6 @@
+"""
+Convenience functions for loading from different sources.
+"""
 
 import json
 import urllib2
@@ -8,8 +11,6 @@ from . import parser
 # USER FUNCTIONS
 #################
 
-# convenience methods for loading from different sources
-
 def from_url(url, format=None):
     """
     Returns the crs object from a string interpreted as a specified format, located at a given url site.

pycrs/parser.py

@@ -1,5 +1,7 @@
+"""
+The main module containing functions for parsing text strings into crs objects.
+"""
 
-# parse from text strings
 # possible use module: https://github.com/rockdoc/grabbag/wiki/CRS-WKT-Parser
 # also note some paramter descriptions: http://www.geoapi.org/3.0/javadoc/org/opengis/referencing/doc-files/WKT.html
 # and see gdal source code: http://gis.stackexchange.com/questions/129764/how-are-esri-wkt-projections-different-from-ogc-wkt-projections
@@ -303,7 +305,8 @@ def _from_wkt(string, wkttype, strict=False):
                 unit = unitclass()
                 unittype = parameters.UnitType(unit)
             else:
-                raise Exception("The specified unit name could not be found")
+                unit = units.Unknown()
+                unittype = parameters.UnitType(unit)
 
             metmult = parameters.MeterMultiplier(value)
             linunit = parameters.Unit(unittype, metmult)
@@ -340,15 +343,16 @@ def _from_wkt(string, wkttype, strict=False):
             if ellipsclass:
                 ellipsdef = ellipsclass()
             else:
-                raise Exception("The specified ellipsoid name could not be found")
+                ellipsdef = ellipsoids.Unknown()
 
             ellipsoid = containers.Ellipsoid(ellipsdef, subsubcontent[1], subsubcontent[2])
 
             ## datum shift
             if wkttype == "ogc":
-                if len(subcontent) >= 3:
-                    subsubheader, subsubcontent = subcontent[2]
-                    datumshift = parameters.DatumShift(subsubcontent)
+                for subsubheader,subsubcontent in subcontent[1:]:
+                    if subsubheader == "TOWGS84":
+                        datumshift = parameters.DatumShift(subsubcontent)
+                        break
                 else:
                     datumshift = None
             elif wkttype == "esri":
@@ -370,7 +374,8 @@ def _from_wkt(string, wkttype, strict=False):
                 unit = unitclass()
                 unittype = parameters.UnitType(unit)
             else:
-                raise Exception("The specified unit name could not be found")
+                unit = units.Unknown()
+                unittype = parameters.UnitType(unit)
             metmult = parameters.MeterMultiplier(value)
             angunit = parameters.AngularUnit(unittype, metmult)
             
@@ -459,12 +464,18 @@ def from_proj4(string, strict=False):
         ellipsclass = ellipsoids.find(ellipsname, "proj4", strict)
         if ellipsclass:
             ellipsdef = ellipsclass
+        elif "+a" in partdict and "+f" in partdict:
+            ellipsdef = ellipsoids.Unknown()
         else:
-            raise Exception("The specified ellipsoid name could not be found")
-
+            raise Exception("The specified ellipsoid name could not be found, and there was no manual specification of the semimajor axis and inverse flattening to use as a substitute.")
 
+    elif "+a" in partdict and "+f" in partdict:
+        # alternatively, it is okay with a missing ellipsoid if +a and +f are specified
+        # TODO: +f seems to never be specified when +ellps is missing, only +a and +b, look into...
+        ellipsdef = ellipsoids.Unknown()
+        
     else:
-        raise Exception("Could not find required +ellps element")
+        raise Exception("Could not find the required +ellps element, nor a manual specification of the +a or +f elements.")
 
     # TO WGS 84 COEFFS
     if "+towgs84" in partdict:

pycrs/utils.py

@@ -1,3 +1,7 @@
+"""
+Misc utility functions related to crs formats and online services. 
+"""
+
 import urllib2
 import re
 

tester.py

@@ -6,6 +6,7 @@ import logging
 
 ###########################
 # Drawing routine for testing
+raw = None
 def render_world(crs, savename):
     import urllib2
     import json
@@ -15,7 +16,9 @@ def render_world(crs, savename):
     import random
 
     # load world borders
-    raw = urllib2.urlopen("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json").read()
+    global raw
+    if not raw:
+        raw = urllib2.urlopen("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json").read()
     rawdict = json.loads(raw)
     data = pygeoj.load(data=rawdict)
 
@@ -215,6 +218,8 @@ def testoutputs(crs):
 
 
 
+#############################################################################
+#############################################################################
 #############################################################################