modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
517 stars 313 forks source link

Reference branch #43

Closed jtwhite79 closed 9 years ago

jtwhite79 commented 9 years ago

Hey guys - I've worked up a feature branch that splits off the spatial reference info from ModelMap into new SpatialReference class. I also created a TemporalReference class for time handling. It took some refactoring of ModelMap to try to separate it from the requirement of a complete flopy.mbase object. I updated the ModelMap notebook and added a new py.test file.

I've already added a plot() and to_shapefile() method to util_2d util_3d transient_2d array so that any array object can be plotted or saved and also added datetime info to mflist so that model input files now have datetime info in the comments. I anticipate this will open up some really cool functionality such as observation pre- and post-processing with pandas, a spatially and temporally referenced binaryfile with plot() methods for both maps and time series, etc...

check it out and let me know what you think. Want to merge soon...

langevin-usgs commented 9 years ago

Jeremy, is the map stuff backward compatible? Will we have to change all of the material for the Python class?

jdhughes-usgs commented 9 years ago

Shouldn't the plot() and to_shapefile() methods be in plot rather than in util? Wouldn't their inclusion in util necessitate a dependency on matplotlib and shapefile?

to_shapefile() should probably be in an export class rather than a method in util_2d, util_3d or transient_2d

Couldn't these methods, if in plot, just take the instance of the util_2d, util_3d or transient_2d and apply the appropriate method (.array)?

Seems a dependency of plot() and to_shapefile() on util_2d, util_3d or transient_2d is better than the other way around.

mbakker7 commented 9 years ago

Hey Guys,

Interesting stuff, Jeremy.

My question is also on backwards compatibility. Is it?

We discussed writing a short (4-page) FloPy paper for Groundwater. I talked to Henk Haitjema (GW Editor), and he thought it was a great idea. I suggest we start a paper on Google Docs and write it from there.

Big question is, of course, what we want to put in. I want to put some examples in, but they should be short enough for people to understand. But I think we also want to talk about all the cool features.

Let me start a doc and share with you, then anybody interested can put in content.

Mark

On Mon, Jul 20, 2015 at 9:31 PM, Hughes, J.D. notifications@github.com wrote:

Shouldn't the plot() and to_shapefile() methods be in plot rather than in util? Wouldn't their inclusion in util necessitate a dependency on matplotlib and shapefile?

to_shapefile() should probably be in an export class rather than a method in util_2d, util_3d or transient_2d

Couldn't these methods, if in plot, just take the instance of the util_2d, util_3d or transient_2d and apply the appropriate method (.array)?

Seems a dependency of plot() and to_shapefile() on util_2d, util_3d or transient_2d is better than the other way around.

— Reply to this email directly or view it on GitHub https://github.com/modflowpy/flopy/issues/43#issuecomment-122999313.

jtwhite79 commented 9 years ago

@langevin-usgs and @mbakker7 - It is backward compatible with one exception - if you don't pass the ax arg to ModelMap, it creates one in the constructor via subplot(111,aspect="equal"). The way it was originally, you could have a ModelMap with self.ax = None, which was causing trouble. But I'm not saying we should use this branch for the class.

@jdhughes-usgs - the plot() and to_shapefile() methods are just thin wrappers around ModelMap. The goal of these is to make flopy more like pandas: a complete set of methods on the core objects for users to explore with. So instead of having to create separate ModelMap instances (plural) with the right args for each one and then call plot_array() or write_grid_shapefile() with the right args for each call, you can just say lpf.hk.plot() since hk already knows what to do - you can get nlay plots of hk (for each layer) with one simple call. Or even better, write the entire recharge sequence to a shapefile with rch.rech.to_shapefile(). And, if you look at the code, you can see that the dependency on ModeMap is only realized if someone calls plot(). I realize this is not very fortran-like - burn!

"haters gonna hate hate hate hate hate"

langevin-usgs commented 9 years ago

shake it off! shake it off!

langevin-usgs commented 9 years ago

Jeremy, what about this. I see that you have put a plot method in the util_array, for example, as shown below in this post. In your plot method, you try to import matplotlib, which is good. That way it doesn't have to be installed. But what if instead of doing it there, all of the plot stuff was moved into the plot utilities. You could still have a very light plot method, similar to what you have, but it would call a function from the plot utilities.

# my proposed approach
def plot(self, axes=None, **kwargs):
    call plot_array(self, axes, kwargs)
    return

If you did it this way, then would could conceivably swap out matplotlib for whatever new plotting libraries become available in the future. If we went this route, then all plotting calls would be made from the plot part of flopy.

# your current implementation
def plot(self, axes=None, **kwargs):
    '''
    How about some doc strings
    '''
    try:
        import matplotlib.pyplot as plt
    except:
        s = 'Could not import matplotlib.  Must install matplotlib ' +\
            ' in order to plot util_3d data.'
        raise Exception(s)
    if axes is not None:
        assert len(axes) == self.shape[0]
        for k in range(self.shape[0]):
            self.util_2ds[k].plot(ax=axes[k], **kwargs)

    else:
        for k in range(self.shape[0]):
            fig = plt.figure()
            ax = plt.subplot(1, 1, 1, aspect='equal')
            self.util_2ds[k].plot(ax=ax, **kwargs)
jdhughes-usgs commented 9 years ago

I have taken a closer look at the reference branch. I made a few changes to make ModelMap behavior more similar to the develop branch. As usual, after taking a look at the functionality I like the way it is implemented. The advantage to the current implementation is the user does not have to know anything about the data to plot model data.

To plot a 2D array the user would enter:

mf.dis.top.plot()

To plot a 3D array the user would enter:

mf.dis.botm.plot()

The implementation is lightweight and intuitive. Furthermore, it accepts matplotlib **kwargs for complete user customization.

A few other useful additions would be 1) creating a LayerFile.plot_data() method and 2) addition of a mflist.plot() method.

A possible implementation of LayerFile.plot_data() method would be:

def plot_data(self, kstpkper=(0, 0), idx=None, totim=0, mflay=None):

which is identical to the LayerFile.get_data method.

The mflist.plot() method could call the ModelMap.plot_bc() method with data for a particular time step in stress_period_data for a package.

I am a convert and suggest pulling it into the develop branch and making a new flopy release with this functionality (and hopefully the proposed functionality) before the San Diego class.

jtwhite79 commented 9 years ago

@langevin-usgs - I like the idea of encapsulating the tedious plotting stuff away from util_array, but I think (and I may be wrong) we may have to come up with separate helpers in plot_util for each of util_2d, util_3d and transient_2d right? Not that its hard to do...I'll work on cleaning it up. I just added the spatial and temporal referencing to LayerFile and have a minimum working code for a plot_data() and to_shapefile() for LayerFile, so now you can do things like

hds = HeadFile(somefile,dis=dis)
hds.plot_data(totim=sometime)
hds.to_shapefile(somefilename,totim=sometime)

This is just a first cut, but it seems to be a useful addition, at least to me.

langevin-usgs commented 9 years ago

Sounds good to me, J. I really like the functionality. Can't wait to use it.

jtwhite79 commented 9 years ago

I added the support for the (soon-to-be) new default colormap in matplotlib as a custom colormap in plotutils. jet is no good: http://betterfigures.org/2015/07/10/a-welcome-development-for-matplotlib/

jdhughes-usgs commented 9 years ago

When it becomes default then it will be the default. I saw we leave it as is for now rather than add it to flopy. As an aside, I think viridis looks like a full diaper,

langevin-usgs commented 9 years ago

Jeremy, I merged reference into develop. I'm getting an error when running py.test, so leaving this comment open for the time being.

jtwhite79 commented 9 years ago

Thanks for merging. I just committed a fix to develop - all tests are passing.