Geografic-Information-System/gui/wxpython/wxplot/profile.py

"""
@package wxplot.profile

@brief Profiling using PyPlot

Classes:
 - profile::ProfileFrame
 - profile::ProfileToolbar

(C) 2011-2014 by the GRASS Development Team

This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.

@author Michael Barton, Arizona State University
"""

import os
import sys
import six
import math
import numpy

import wx

import wx.lib.plot as plot
import grass.script as grass
from wxplot.base import BasePlotFrame, PlotIcons
from gui_core.toolbars import BaseToolbar, BaseIcons
from gui_core.wrap import StockCursor
from wxplot.dialogs import ProfileRasterDialog, PlotStatsFrame
from core.gcmd import RunCommand, GWarning, GError, GMessage


class ProfileFrame(BasePlotFrame):
    """Mainframe for displaying profile of one or more raster maps. Uses wx.lib.plot.
    """

    def __init__(self, parent, giface, controller, units, size=wx.Size(700, 400),
                 rasterList=None, **kwargs):
        BasePlotFrame.__init__(self, parent=parent, giface=giface, size=size, **kwargs)

        self.controller = controller
        self.controller.transectChanged.connect(self.SetTransect)
        self.transect = []
        self.toolbar = ProfileToolbar(parent=self)
        # workaround for http://trac.wxwidgets.org/ticket/13888
        if sys.platform != 'darwin':
            self.SetToolBar(self.toolbar)
        self.SetTitle(_("GRASS Profile Analysis Tool"))
        # in case of degrees, we want to use meters
        self._units = units if 'degree' not in units else 'meters'

        #
        # Init variables
        #
        if rasterList is None:
            self.rasterList = []
        else:
            self.rasterList = rasterList
        self.plottype = 'profile'
        self.coordstr = ''              # string of coordinates for r.profile
        self.seglist = []               # segment endpoint list
        self.ppoints = ''               # segment endpoints data
        self.transect_length = 0.0      # total transect length
        self.ptitle = _('Profile of')   # title of window
        self.colorList = [
            "blue",
            "red",
            "green",
            "yellow",
            "magenta",
            "cyan",
            "aqua",
            "black",
            "grey",
            "orange",
            "brown",
            "purple",
            "violet",
            "indigo"]

        self._initOpts()

        if len(
                self.rasterList) > 0:  # set raster name(s) from layer manager if a map is selected
            self.raster = self.InitRasterOpts(self.rasterList, self.plottype)
        else:
            self.raster = {}

        # determine units (axis labels)
        # maybe, we should not accept these invalid units
        # but ok, trying to handle it here
        if self._units is not None and self._units != '':
            self.xlabel = _('Distance (%s)') % self._units
        else:
            self.xlabel = _("Distance along transect")
        self.ylabel = _("Cell values")

        # Bind events
        self.Bind(wx.EVT_CLOSE, self.OnCloseWindow)

        self.SetGraphStyle()

    def _initOpts(self):
        """Initialize plot options
        """
        self.InitPlotOpts('profile')

    def SetTransect(self, coords):
        self.transect = coords
        if coords:
            self.OnCreateProfile(None)
        else:
            self.OnErase(None)

    def OnDrawTransect(self, event):
        """Draws transect to profile in map display
        """
        if self.controller.IsActive():
            self.controller.Stop()
        self.controller.Start()

        self.parent.SetFocus()
        self.parent.Raise()

    def OnSelectRaster(self, event):
        """Select raster map(s) to profile
        """
        dlg = ProfileRasterDialog(parent=self)
        dlg.CenterOnParent()
        if dlg.ShowModal() == wx.ID_OK:
            self.rasterList = dlg.rasterList
            self.raster = self.InitRasterOpts(self.rasterList, self.plottype)

            # plot profile
            if len(self.transect) > 0 and len(self.rasterList) > 0:
                self.OnCreateProfile(event=None)

        dlg.Destroy()

    def SetupProfile(self):
        """Create coordinate string for profiling. Create segment
           list for transect segment markers.
        """
        # create list of coordinate points for r.profile
        dist = 0
        cumdist = 0
        self.coordstr = ''
        lasteast = lastnorth = None

        region = grass.region()
        insideRegion = True
        if len(self.transect) > 0:
            for point in self.transect:
                if not (region['w'] <= point[0] <= region['e']
                        and region['s'] <= point[1] <= region['n']):
                    insideRegion = False
                # build string of coordinate points for r.profile
                if self.coordstr == '':
                    self.coordstr = '%f,%f' % (point[0], point[1])
                else:
                    self.coordstr = '%s,%f,%f' % (
                        self.coordstr, point[0], point[1])

        if not insideRegion:
            GWarning(
                message=_(
                    "Not all points of profile lie inside computational region."),
                parent=self)

        if len(self.rasterList) == 0:
            return

        # title of window
        self.ptitle = _('Profile of')

        # create list of coordinates for transect segment markers
        if len(self.transect) > 0:
            self.seglist = []
            for point in self.transect:
                # get value of raster cell at coordinate point
                ret = RunCommand('r.what',
                                 parent=self,
                                 read=True,
                                 map=self.rasterList[0],
                                 coordinates='%f,%f' % (point[0], point[1]))

                val = ret.splitlines()[0].split('|')[3]
                if val is None or val == '*':
                    continue
                val = float(val)

                # calculate distance between coordinate points
                if lasteast and lastnorth:
                    dist = math.sqrt(
                        math.pow((lasteast - point[0]),
                                 2) + math.pow((lastnorth - point[1]),
                                               2))
                cumdist += dist

                # store total transect length
                self.transect_length = cumdist

                # build a list of distance,value pairs for each segment of
                # transect
                self.seglist.append((cumdist, val))
                lasteast = point[0]
                lastnorth = point[1]

            # delete extra first segment point
            try:
                self.seglist.pop(0)
            except:
                pass

        #
        # create datalist of dist/value pairs and y labels for each raster map
        #
        self.ylabel = ''
        i = 0

        for r in six.iterkeys(self.raster):
            self.raster[r]['datalist'] = []
            datalist = self.CreateDatalist(r, self.coordstr)
            if len(datalist) > 0:
                self.raster[r]['datalist'] = datalist

                # update ylabel to match units if they exist
                if self.raster[r]['units'] != '':
                    self.ylabel += '%s (%d),' % (self.raster[r]['units'], i)
                i += 1

                # update title
                self.ptitle += ' %s ,' % r.split('@')[0]

        self.ptitle = self.ptitle.rstrip(',')

        if self.ylabel == '':
            self.ylabel = _('Raster values')
        else:
            self.ylabel = self.ylabel.rstrip(',')

    def CreateDatalist(self, raster, coords):
        """Build a list of distance, value pairs for points along transect using r.profile
        """
        datalist = []

        # keep total number of transect points to 500 or less to avoid
        # freezing with large, high resolution maps
        region = grass.region()
        curr_res = min(float(region['nsres']), float(region['ewres']))
        transect_rec = 0
        if self.transect_length / curr_res > 500:
            transect_res = self.transect_length / 500
        else:
            transect_res = curr_res

        ret = RunCommand("r.profile",
                         parent=self,
                         input=raster,
                         coordinates=coords,
                         resolution=transect_res,
                         null="nan",
                         quiet=True,
                         read=True)

        if not ret:
            return []

        for line in ret.splitlines():
            dist, elev = line.strip().split(' ')
            if dist is None or dist == '' or dist == 'nan' or \
                    elev is None or elev == '' or elev == 'nan':
                continue
            dist = float(dist)
            elev = float(elev)
            datalist.append((dist, elev))

        return datalist

    def OnCreateProfile(self, event):
        """Main routine for creating a profile. Uses r.profile to
        create a list of distance,cell value pairs. This is passed to
        plot to create a line graph of the profile. If the profile
        transect is in multiple segments, these are drawn as
        points. Profile transect is drawn, using methods in mapdisp.py
        """

        if len(self.transect) == 0 or len(self.rasterList) == 0:
            dlg = wx.MessageDialog(
                parent=self,
                message=_('You must draw a transect to profile in the map display window.'),
                caption=_('Nothing to profile'),
                style=wx.OK | wx.ICON_INFORMATION | wx.CENTRE)
            dlg.ShowModal()
            dlg.Destroy()
            return

        self.SetCursor(StockCursor(wx.CURSOR_ARROW))

        self.SetupProfile()
        p = self.CreatePlotList()
        self.DrawPlot(p)

    def CreatePlotList(self):
        """Create a plot data list from transect datalist and
            transect segment endpoint coordinates.
        """
        # graph the distance, value pairs for the transect
        self.plotlist = []

        # Add segment marker points to plot data list
        if len(self.seglist) > 0:
            self.ppoints = plot.PolyMarker(
                self.seglist,
                legend=' ' + self.properties['marker']['legend'],
                colour=wx.Colour(
                    self.properties['marker']['color'][0],
                    self.properties['marker']['color'][1],
                    self.properties['marker']['color'][2],
                    255),
                size=self.properties['marker']['size'],
                fillstyle=self.ptfilldict[
                    self.properties['marker']['fill']],
                marker=self.properties['marker']['type'])
            self.plotlist.append(self.ppoints)

        # Add profile distance/elevation pairs to plot data list for each
        # raster profiled
        for r in self.rasterList:
            col = wx.Colour(self.raster[r]['pcolor'][0],
                            self.raster[r]['pcolor'][1],
                            self.raster[r]['pcolor'][2],
                            255)
            self.raster[r]['pline'] = plot.PolyLine(
                self.raster[r]['datalist'],
                colour=col, width=self.raster[r]['pwidth'],
                style=self.linestyledict[self.raster[r]['pstyle']],
                legend=self.raster[r]['plegend'])

            self.plotlist.append(self.raster[r]['pline'])

        if len(self.plotlist) > 0:
            return self.plotlist
        else:
            return None

    def Update(self):
        """Update profile after changing options
        """
        self.SetGraphStyle()
        p = self.CreatePlotList()
        self.DrawPlot(p)

    def SaveProfileToFile(self, event):
        """Save r.profile data to a csv file
        """
        dlg = wx.FileDialog(
            parent=self,
            message=_(
                "Choose prefix for file(s) where to save profile values..."),
            defaultDir=os.getcwd(),
            wildcard=_("Comma separated value (*.csv)|*.csv"),
            style=wx.FD_SAVE)
        pfile = []
        if dlg.ShowModal() == wx.ID_OK:
            path = dlg.GetPath()
            for r in self.rasterList:
                pfile.append(path + '_' + str(r.replace('@', '_')) + '.csv')
                if os.path.exists(pfile[-1]):
                    dlgOv = wx.MessageDialog(
                        self,
                        message=_(
                            "File <%s> already exists. "
                            "Do you want to overwrite this file?") % pfile
                        [-1],
                        caption=_("Overwrite file?"),
                        style=wx.YES_NO | wx.YES_DEFAULT | wx.ICON_QUESTION)
                    if dlgOv.ShowModal() != wx.ID_YES:
                        pfile.pop()
                        dlgOv.Destroy()
                        continue

                try:
                    fd = open(pfile[-1], "w")
                except IOError as e:
                    GError(parent=self,
                           message=_("Unable to open file <%s> for writing.\n"
                                     "Reason: %s") % (pfile[-1], e))
                    dlg.Destroy()
                    return

                for datapair in self.raster[r]['datalist']:
                    fd.write(
                        '%.6f,%.6f\n' %
                        (float(
                            datapair[0]), float(
                            datapair[1])))

                fd.close()

        dlg.Destroy()
        if pfile:
            message = _("%d files created:\n%s") % (
                len(pfile), '\n'.join(pfile))
        else:
            message = _("No files generated.")

        GMessage(parent=self, message=message)

    def OnStats(self, event):
        """Displays regression information in messagebox
        """
        message = []
        title = _('Statistics for Profile(s)')

        for r in six.iterkeys(self.raster):
            try:
                rast = r.split('@')[0]
                statstr = 'Profile of %s\n\n' % rast

                iterable = (i[1] for i in self.raster[r]['datalist'])
                a = numpy.fromiter(iterable, numpy.float)

                statstr += 'n: %f\n' % a.size
                statstr += 'minimum: %f\n' % numpy.amin(a)
                statstr += 'maximum: %f\n' % numpy.amax(a)
                statstr += 'range: %f\n' % numpy.ptp(a)
                statstr += 'mean: %f\n' % numpy.mean(a)
                statstr += 'standard deviation: %f\n' % numpy.std(a)
                statstr += 'variance: %f\n' % numpy.var(a)
                cv = numpy.std(a) / numpy.mean(a)
                statstr += 'coefficient of variation: %f\n' % cv
                statstr += 'sum: %f\n' % numpy.sum(a)
                statstr += 'median: %f\n' % numpy.median(a)
                statstr += 'distance along transect: %f\n\n' % self.transect_length
                message.append(statstr)
            except:
                pass

        stats = PlotStatsFrame(self, id=wx.ID_ANY, message=message,
                               title=title)

        if stats.Show() == wx.ID_CLOSE:
            stats.Destroy()

    def OnCloseWindow(self, event):
        if self.controller.IsActive():
            self.controller.Stop()
        self.Destroy()


class ProfileToolbar(BaseToolbar):
    """Toolbar for profiling raster map
    """

    def __init__(self, parent):
        BaseToolbar.__init__(self, parent)

        # workaround for http://trac.wxwidgets.org/ticket/13888
        if sys.platform == 'darwin':
            parent.SetToolBar(self)

        self.InitToolbar(self._toolbarData())

        # realize the toolbar
        self.Realize()

    def _toolbarData(self):
        """Toolbar data"""
        return self._getToolbarData((('addraster', BaseIcons["addRast"],
                                      self.parent.OnSelectRaster),
                                     ('transect', PlotIcons["transect"],
                                      self.parent.OnDrawTransect),
                                     (None, ),
                                     ('draw', PlotIcons["draw"],
                                      self.parent.OnCreateProfile),
                                     ('erase', BaseIcons["erase"],
                                      self.parent.OnErase),
                                     ('drag', BaseIcons['pan'],
                                      self.parent.OnDrag),
                                     ('zoom', BaseIcons['zoomIn'],
                                      self.parent.OnZoom),
                                     ('unzoom', BaseIcons['zoomBack'],
                                      self.parent.OnRedraw),
                                     (None, ),
                                     ('statistics', PlotIcons['statistics'],
                                      self.parent.OnStats),
                                     ('datasave', PlotIcons["save"],
                                      self.parent.SaveProfileToFile),
                                     ('image', BaseIcons["saveFile"],
                                      self.parent.SaveToFile),
                                     ('print', BaseIcons["print"],
                                      self.parent.PrintMenu),
                                     (None, ),
                                     ('settings', PlotIcons["options"],
                                      self.parent.PlotOptionsMenu),
                                     ('quit', PlotIcons["quit"],
                                      self.parent.OnQuit),
                                     ))