Source code for geophpy.geoposset

# -*- coding: utf-8 -*-
'''
    GeophPy.geoposset
    -----------------

    Geographics Positioning Set management.

    :copyright: Copyright 2014 Lionel Darras, Philippe Marty and contributors, see AUTHORS.
    :license: GNU GPL v3.

'''
from __future__ import absolute_import
import shapefile   # pyshp for reading and writing shapefiles
import geophpy.geopositioning.kml as kml
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import glob                         # for managing severals files thanks to "*." extension
import utm                          # to convert utm and wgs84 coordinates
import os

# definitions
POINT_PARAMS = 'bo'
MULTIPOINT_PARAMS = 'bo--'
LINE_PARAMS = 'r'
SURFACE_PARAMS = 'g'

refsystem_list = ['UTM', 'WGS84']           # list of reference systems availables
type_list = ['point', 'shapefile']          # list of element types availables
utm_minletter = 'E'
utm_minnumber = 1
utm_maxletter = 'X'
utm_maxnumber = 60


#---------------------------------------------------------------------------#
# Geographics Positioning Set Object                                        #
#---------------------------------------------------------------------------#
[docs]class GeoPosSet(object): ''' Builds a geographic positioning set to display or use to georeferencing data set. ''' refsystem = None # Geographic positioning Reference System (None, 'UTM', 'WGS84', ...) type = None # Type of elements contained in the array ('point', 'shapefile', ...) list = [] # List of elements (points, shapefiles, ...) def __init__(self): pass @classmethod
[docs] def from_file(cls, refsystem, type, filenameslist): '''Builds a geographic positioning set from a list of files. Parameters : :refsystem: geographic positioning reference system, None, 'UTM', 'WGS84', ... :type: type of elements contained in the files, 'point', 'shapefile', ... :filenameslist: list of files to open, ['shapefile1', 'shapefile2', ...] or, ['shapefile*'] to open all files with a filename beginning by "shapefile" (without extensions). Returns : :success: true if geoposset built, fals if not. :geoposset: geographics positioning set built with file(s). Example : success, geoposset = GeoPosSet.from_file('UTM', 'shapefile', "shapefile1") ''' geoposset = cls() geoposset.type = type geoposset.refsystem = refsystem geoposset.list = [] success = True # by default if (type == 'shapefile') : for filenames in filenameslist : # for each file in the list try: filenamewithoutextension, extension = os.path.splitext(filenames) filenamesextended = glob.glob(filenamewithoutextension+".shp") # extension if the filenames if '*' is contained for filename in filenamesextended : geoposset.list.append(shapefile.Reader(filename)) # adds shapefile in the list except IOError: print("File Read Error (",filename,")") success = False break # exits the loop else : success = False return success, geoposset
[docs] def points_getlist(self): ''' get list of points Returns: list of points ''' point_num = 0 # initialisation of the point number points_list = [] # initialisation of the points list if (self.type == 'shapefile'): firstpoint = True for sf in self.list : # for each shape file in the list shapes = sf.shapes() for shape in shapes : # for each shape in the file # if point type : if ((shape.shapeType == shapefile.POINT) or (shape.shapeType == shapefile.POINTZ) or (shape.shapeType == shapefile.POINTM)) : points_list.append([point_num, shape.points[0][0], shape.points[0][1]]) point_num = point_num + 1 else : points_list = [] return points_list
[docs] def plot(self, picturefilename=None, dpi=None, transparent=False, i_xmin=None, i_xmax=None, i_ymin=None, i_ymax=None): ''' display map of points, lines, and surfaces described in geographic positioning list Parameters : :picturefilename: name of the file to save the picture. If None, no picture saved in a file :dpi: :transparent: if True, picture display points not plotted as transparents :i_xmin: x minimal value to display, None by default :i_xmax: x maximal value to display, None by default :i_ymin: y minimal value to display, None by default :i_ymax: y maximal value to display, None by default Returns : :success: True if no error :fig: figure object ''' fig = plt.figure() # creation of the empty figure point_num = 0 # initialisation of the point number points_list = [["Num", "Lon", "Lat"]] # initialisation of the points list firstpoint = True success = True if (self.type == 'shapefile'): for sf in self.list : # for each shape file in the list shapes = sf.shapes() for shape in shapes : # for each shape in the file # if point type : if ((shape.shapeType == shapefile.POINT) or (shape.shapeType == shapefile.POINTZ) or (shape.shapeType == shapefile.POINTM)) : plot_params = MULTIPOINT_PARAMS points_list.append([point_num, shape.points[0][0], shape.points[0][1]]) plt.text(shape.points[0][0], shape.points[0][1],str(point_num)) point_num = point_num + 1 # if multipoint type : elif ((shape.shapeType == shapefile.MULTIPOINT) or (shape.shapeType == shapefile.MULTIPOINTZ) or (shape.shapeType == shapefile.MULTIPOINTM)) : plot_params = MULTIPOINT_PARAMS # if line type : elif ((shape.shapeType == shapefile.POLYLINE) or (shape.shapeType == shapefile.POLYLINEZ) or (shape.shapeType == shapefile.POLYLINEM)) : plot_params = LINE_PARAMS # if surface type : else : plot_params = SURFACE_PARAMS points = np.array(shape.points) if (firstpoint == True) : # if first point to display, initialises the x and y limits values xmin = min(points.T[0]) xmax = max(points.T[0]) ymin = min(points.T[1]) ymax = max(points.T[1]) firstpoint = False else : # else compares with the x and y limits established if (min(points.T[0]) < xmin) : xmin = min(points.T[0]) if (max(points.T[0]) > xmax) : xmax = max(points.T[0]) if (min(points.T[1]) < ymin) : ymin = min(points.T[1]) if (max(points.T[1]) > ymax) : ymax = max(points.T[1]) # displays the point, the line, or the surface plt.plot(points.T[0], points.T[1], plot_params) else: success = False if (success == True): # for each x or y limit not configured in input, initialisation of the value with the limits of points displayed nolimits = True if (i_xmin == None) : i_xmin = xmin else : nolimits = False if (i_xmax == None) : i_xmax = xmax else : nolimits = False if (i_ymin == None) : i_ymin = ymin else : nolimits = False if (i_ymax == None) : i_ymax = ymax else : nolimits = False if (nolimits == False) : # sets the axis limits xmin, xmax, ymin, ymax = plt.axis([i_xmin, i_xmax, i_ymin, i_ymax]) else : # gets the axis limits xmin, xmax, ymin, ymax = plt.axis() # to have the same scale in X and Y axis dy = ymax - ymin dx = xmax - xmin if (dy > dx): xmax = xmin + dy elif (dx > dy): ymax = ymin + dx plt.axis([xmin, xmax, ymin, ymax]) ax = plt.gca() ax.ticklabel_format(useOffset=False) if (picturefilename != None) : # if picturefilename in input, saves in a picture file plt.savefig(picturefilename, dpi=dpi, transparent=transparent) return success, fig
[docs] def to_kml(self, kmlfilename): """ Builds a kml file with points, lines, and surfaces described in a list of geographic positioning set Parameters : :kmlfilename: Name of the kml file to save Returns : :success: True if operation done """ success = True # by default if (self.type == 'shapefile'): kml.shapefile_to_kmlfile(self.list, kmlfilename) else: success = False return success #===========================================================================# #---------------------------------------------------------------------------# # General Geographic Positioning Set Environments functions # #---------------------------------------------------------------------------#
[docs]def refsys_getlist(): """ Returns list of reference systems availables, 'UTM', 'WGS84', ... """ return refsystem_list
[docs]def type_getlist(): """ Returns list of geographic positioning type, 'point', 'shapefile', ... """ return type_list
[docs]def utm_to_wgs84(easting, northing, zonenumber, zoneletter): """ Converts UTM to WGS84 coordinates Parameters : :easting: east UTM coordinate :northing: north UTM coordinate :zonenumber: zone number :zone letter: zone letter Returns : :latitude: WGS84 latitude coordinate :longitude: WGS84 longitude coordinate """ return utm.to_latlon(easting, northing, zonenumber, zoneletter)
[docs]def wgs84_to_utm(latitude, longitude): """ Converts WGS84 to UTM coordinates Parameters : :latitude: WGS84 latitude coordinate :longitude: WGS84 longitude coordinate Returns : :easting: east UTM coordinate :northing: north UTM coordinate :zonenumber: zone number :zone letter: zone letter """ return utm.from_latlon(latitude, longitude)
[docs]def utm_getzonelimits(): """ Gets the min and max numbers and letters of the UTM coordinates system Returns: :min_number: minimal number of the UTM zone, 1. :min_letter: minimal letter of the UTM zone, E. :max_numbre: maximal number of the UTM zone, 60. :max letter: maximal letter of the UTM zone, X. """ return utm_minnumber, utm_minletter, utm_maxnumber, utm_maxletter