<jupyter_start><jupyter_text>FloPy  `SpatialReference` demoA short demonstration of functionality in the `SpatialReference` class for locating the model in a "real world" coordinate reference system<jupyter_code>import sys
import os
import shutil
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# run installed version of flopy or add local path
try:
    import flopy
except:
    fpth = os.path.abspath(os.path.join('..', '..'))
    sys.path.append(fpth)
    import flopy
    
from flopy.utils.reference import SpatialReference
import flopy.utils.binaryfile as bf

print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('flopy version: {}'.format(flopy.__version__))
outpath = 'data'
if os.path.isdir(outpath):
    shutil.rmtree(outpath)
os.mkdir(outpath)<jupyter_output><empty_output><jupyter_text>Description`SpatialReference` is a stand-alone option that describes the location of the model grid in a real-world coordinate reference system. It can be created with some basic information about the model's location.<jupyter_code>nrow, ncol = 10, 10
xll, yll = 617822.3, 5114820.7 # origin of the model (lower left corner)
dxdy = 250 # grid spacing (in model units)
rot = 29 # rotation (positive counterclockwise)

# epsg code specifying coordinate reference system
# (https://www.epsg-registry.org/)
# in this case, UTM zone 16 N, NAD83
# https://spatialreference.org/ref/epsg/nad83-utm-zone-16n/
model_epsg = 26915

# alternatively, a proj4 string can be supplied
model_proj4 = 'https://spatialreference.org/ref/epsg/nad83-utm-zone-16n/proj4/'

# row and column spacings
# (note that delc is column spacings along a row; delr the row spacings along a column)
delc = np.ones(nrow, dtype=float) * dxdy
delr = np.ones(ncol, dtype=float) * dxdy

sr = SpatialReference(delr=delr, delc=delc, xll=xll, yll=yll, rotation=rot, epsg=model_epsg)<jupyter_output><empty_output><jupyter_text>`SpatialReference` makes many calculations internally, so you don't have tocalling the sr object prints the important parameters. Note that the upper right corner of the model, which is often used to specify the origin, has been calculated. A proj4 string has also been fetched from  using the supplied epsg code.<jupyter_code>sr<jupyter_output><empty_output><jupyter_text>Note the length parameters**units**: Length units for real world coordinate system (typically meters; feet are also supported). This attribute is inferred from the epsg code or proj4 string, but can also be supplied.  **lenuni**: MODFLOW length unit ([see documentation]()) (default is 2, meters).  **length_multiplier**: multiplier for scaling grid from MODFLOW units to real world crs units. This parameter is inferred from the above units, or can be supplied directly.   `SpatialReference` parameters can be updated dynamically<jupyter_code># switch the MODFLOW units to feet
sr.lenuni = 1
sr<jupyter_output><empty_output><jupyter_text>Note that the length multiplier and xul, yul changedthe lower left corner was maintained because it was specified. If xul, yul had been specified instead, that would have been preserved.<jupyter_code>sr.xll, sr.yll<jupyter_output><empty_output><jupyter_text>information about the grid coordinates cell centroids<jupyter_code>sr.xcentergrid, sr.ycentergrid<jupyter_output><empty_output><jupyter_text>cell vertices<jupyter_code>sr.vertices<jupyter_output><empty_output><jupyter_text>Transformations model coordinates to real-world coordinatese.g., for working with MODPATH output<jupyter_code>sr.transform(np.array([0, 10.]), np.array([0, 2.5]))<jupyter_output><empty_output><jupyter_text>real-word coordinates to model coordinatese.g., for specifying observations for Hydmod<jupyter_code>sr.transform(617824.59, 5114822.84, inverse=True)<jupyter_output><empty_output><jupyter_text>model row, column for real-world coordinate<jupyter_code>sr.get_ij(671900., 618000.)<jupyter_output><empty_output><jupyter_text>Grid bounds<jupyter_code>sr.bounds<jupyter_output><empty_output><jupyter_text>Shapefile of grid<jupyter_code>sr.write_shapefile(os.path.join(outpath, 'grid.shp'))<jupyter_output>wrote data/grid.shp<jupyter_text>grid spec file(for PEST)<jupyter_code>sr.write_gridSpec(os.path.join(outpath, 'grid.spc'))<jupyter_output><empty_output><jupyter_text>Loading `SpatialReference``SpatialReference` is included in flopy models at the model level;  with `model.write()` the basic parameters are written to the comment header of the Name file.<jupyter_code># load an existing model
model_ws = os.path.join("..", "data", "freyberg_multilayer_transient")
with open(model_ws + '/freyberg.nam') as input:
    print(input.readlines()[0:2])<jupyter_output>['# Name file for mf2005; generated by Flopy.\n', '#xul:619653; yul:3353277; rotation:15; proj4_str:+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs; start_datetime:1/1/2015\n']<jupyter_text>on load`SpatialReference` info is read from  1) **usgs.model.reference**, if it exists (https://water.usgs.gov/ogw/policy/gw-model/modelers-setup.html)  2) otherwise, the Name file  if no spatial reference information is found, a default sr object is created<jupyter_code>ml = flopy.modflow.Modflow.load("freyberg.nam", model_ws=model_ws, verbose=False,
                               check=False, exe_name="mfnwt")
ml.sr<jupyter_output><empty_output><jupyter_text>Adding/modifying `SpatialReference` for a model<jupyter_code>ml.sr = SpatialReference(delr=ml.dis.delr, delc=ml.dis.delc, 
                         xul=150000, yul=30000, rotation=25, epsg=26715)
ml.sr<jupyter_output><empty_output><jupyter_text>individual parameters can also be updated<jupyter_code>ml.sr.proj4_str = "+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
ml.sr<jupyter_output><empty_output><jupyter_text>note that the epsg code was cleared when the proj4 string was updated<jupyter_code>ml.sr.epsg<jupyter_output><empty_output><jupyter_text>interpolation between cells(e.g. for head values at point locations) read in some head results<jupyter_code>hdsobj = bf.HeadFile(model_ws + '/freyberg.hds')
hds = hdsobj.get_data()
hds[hds < -999] = np.nan
plt.imshow(hds[0])
plt.colorbar();
ml.sr.bounds<jupyter_output><empty_output><jupyter_text>get interpolated head at some pointsuses `scipy.interpolate.gridddata` to interpolate values between cell centers<jupyter_code>ml.sr.interpolate(hds[0], ([157000, 160000], [30000, 27000]))
ml.sr.xll = 0
ml.sr.yll = 0
ml.sr.rotation = 0<jupyter_output><empty_output><jupyter_text>Rasterizing featuresrequires the `rasterio` and `fiona` modules, available via **conda-forge** (mac or windows), or the **Unofficial Python Binaries** website (windows) create a geoJSON-style geometry(this could be easily read in from a shapefile<jupyter_code>try:
    import fiona
except:
    fiona = None
    print('fiona is not available')
try:
    from rasterio.features import rasterize
    from rasterio import Affine
except:
    rasterize = None
    Affine = None
    print('rasterio is not available')<jupyter_output><empty_output><jupyter_text>read in a shapefile<jupyter_code>if fiona is not None:
    shpname = '../data/freyberg/gis/bedrock_outcrop_hole.shp'
    with fiona.open(shpname) as src:
        records = [r for r in src]<jupyter_output><empty_output><jupyter_text>list the feature geometries and attributes<jupyter_code>if fiona is not None:
    geoms = [r['geometry'] for r in records]
    attr = [r['properties'] for r in records]
    attr<jupyter_output><empty_output><jupyter_text>select an attribute to map onto the model grid<jupyter_code>if fiona is not None:
    geoms = [(g, attr[i]['OBJECTID']) for i, g in enumerate(geoms)]<jupyter_output><empty_output><jupyter_text>make an `rasterio.Affine` object using info from the sr<jupyter_code>if rasterize is not None:
    dx = ml.dis.delr.array[0]
    dy = ml.dis.delc.array[0]
    trans = Affine(dx, ml.sr.rotation, ml.sr.xul, 
                   ml.sr.rotation, -dy, ml.sr.yul)
    trans<jupyter_output><empty_output><jupyter_text>rasterize the features, applying the attribute values to cells intersecting each oneIn this case, the cells with zeros are those that don't intersect a feature<jupyter_code>if rasterize is not None:
    r = rasterize(geoms, out_shape=(ml.nrow, ml.ncol), transform=trans)

    fig, ax = plt.subplots()
    qm = ml.sr.plot_array(r, ax=ax)
    ax.set_aspect(1)
    plt.colorbar(qm)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10699539931579569, 'fraction_numerical': 0.020997994573551965, 'mean_word_length': 5.809638554216868, 'pattern_counts': {'":': 0, '<': 104, '<?xml version=': 0, '>': 103, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '87384', 'license': 'CC0-1.0', 'n_tokens_mistral': 3072, 'n_tokens_neox': 2884, 'n_words': 763, 'path': 'examples/Notebooks/flopy3_SpatialReference_demo.ipynb', 'repo_name': 'aleaf/flopy', 'size': '8477'}