Package utils

Shapefile Creator

Description:

Create shape files for Phase Tensor Ellipses, Tipper Real/Imag. export the phase tensor map and tippers into jpeg/png images

CreationDate: 2017-03-06 Developer: fei.zhang@ga.gov.au

Revision History:

LastUpdate: 10/11/2017 FZ fix bugs after the big merge LastUpdate: 20/11/2017 change from freq to period filenames, allow to specify a period LastUpdate: 30/10/2018 combine ellipses and tippers together, refactorings

brenainn.moushall@ga.gov.au 27-03-2020 17:33:23 AEDT:

Fix outfile/directory issue (see commit messages)

class mtpy.utils.shapefiles_creator.ShapefilesCreator(edifile_list, outdir, epsg_code=4326)

Extend the EdiCollection parent class, create phase tensor and tipper shapefiles for a list of edifiles

Parameters

• edifile_list – [path2edi,…]

• outdir – path2output dir, where the shp file will be written.

• {'init' (orig_crs =) – ‘epsg:4283’} # GDA94

Methods

calculate_aver_impedance(dest_dir[, ...])	calculate the average impedance tensor Z (related to apparent resistivity) of all edi (MT-stations) for each period.
create_measurement_csv(dest_dir[, ...])	create csv file from the data of EDI files: IMPEDANCE, APPARENT RESISTIVITIES AND PHASES see also utils/shapefiles_creator.py
create_mt_station_gdf([outshpfile])	create station location geopandas dataframe, and output to shape file
create_penetration_depth_csv(dest_dir[, ...])	create penetration depth csv file for each frequency corresponding to the given input 1.0/period_list.
create_phase_tensor_csv(dest_dir[, ...])	create phase tensor ellipse and tipper properties.
create_phase_tensor_csv_with_image(dest_dir)	Using PlotPhaseTensorMaps class to generate csv file of phase tensor attributes, etc.
create_phase_tensor_shp(period[, ellipsize, ...])	create phase tensor ellipses shape file correspond to a MT period :return: (geopdf_obj, path_to_shapefile)
create_tipper_imag_shp(period[, ...])	create imagery tipper lines shapefile from a csv file The shapefile consists of lines without arrow.
create_tipper_real_shp(period[, ...])	create real tipper lines shapefile from a csv file The shapefile consists of lines without arrow.
display_on_basemap()	display MT stations which are in stored in geopandas dataframe in a base map.
display_on_image()	display/overlay the MT properties on a background geo-referenced map image
export_edi_files(dest_dir[, period_list, ...])	export edi files. :param dest_dir: output directory :param period_list: list of periods; default=None, in which data for all available frequencies are output :param interpolate: Boolean to indicate whether to interpolate data onto given period_list; otherwise a period_list is obtained from get_periods_by_stats() :param file_name: output file name :param period_buffer: buffer so that interpolation doesn't stretch too far over periods. Provide a float or integer factor, greater than which interpolation will not stretch. e.g. 1.5 means only interpolate to a maximum of 1.5 times each side of each frequency value.
get_bounding_box([epsgcode])	compute bounding box
get_min_max_distance()	get the min and max distance between all possible pairs of stations.
get_period_occurance(aper)	For a given aperiod, compute its occurance frequencies among the stations/edi :param aper: a float value of the period :return:
get_periods_by_stats([percentage])	check the presence of each period in all edi files, keep a list of periods which are at least percentage present :return: a list of periods which are present in at least percentage edi files
get_phase_tensor_tippers(period[, interpolate])	For a given MT period (s) value, compute the phase tensor and tippers etc.
get_station_utmzones_stats()	A simple method to find what UTM zones these (edi files) MT stations belong to are they in a single UTM zone, which corresponds to a unique EPSG code? or do they belong to multiple UTM zones?
get_stations_distances_stats()	get the min max statistics of the distances between stations.
plot_stations([savefile, showfig])	Visualise the geopandas df of MT stations
select_periods([show, period_list, percentage])	Use edi_collection to analyse the whole set of EDI files
show_obj([dest_dir])	test call object's methods and show it's properties

create_phase_tensor_shp(period, ellipsize=None, target_epsg_code=4283, export_fig=False)

create phase tensor ellipses shape file correspond to a MT period :return: (geopdf_obj, path_to_shapefile)

create_tipper_imag_shp(period, line_length=None, target_epsg_code=4283, export_fig=False)

create imagery tipper lines shapefile from a csv file The shapefile consists of lines without arrow. User can use GIS software such as ArcGIS to display and add an arrow at each line’s end line_length is how long will be the line, auto-calculatable :return:(geopdf_obj, path_to_shapefile)

create_tipper_real_shp(period, line_length=None, target_epsg_code=4283, export_fig=False)

create real tipper lines shapefile from a csv file The shapefile consists of lines without arrow. User can use GIS software such as ArcGIS to display and add an arrow at each line’s end line_length is how long will be the line, auto-calculatable

mtpy.utils.shapefiles_creator.create_ellipse_shp_from_csv(csvfile, esize=0.03, target_epsg_code=4283)

create phase tensor ellipse geometry from a csv file. This function needs csv file as its input. :param csvfile: a csvfile with full path :param esize: ellipse size, defaut 0.03 is about 3KM in the max ellipse rad :return: a geopandas dataframe

mtpy.utils.shapefiles_creator.create_tensor_tipper_shapefiles(edi_dir, out_dir, periods, pt_base_size=None, pt_phi_max=None, src_epsg=4326, dst_epsg=4326)

Interface for creating and saving phase tensor and tipper shapefiles.

Parameters

edi_dirstr

Path to directory containing .edi data files.

out_dirstr

Path to directory to save resulint shapefiles.

src_epsgint

EPSG code of the EDI data CRS. Defaults 4326 (WGS84).

dst_epsgint

EPSG code of the output (i.e. same CRS as the geotiff you will be displaying on). Defaults 4326 (WGS84).

period_indiciesfloat or list of float. Defaults to 0.0.

List of periods in seconds to create shapefiles for. The nearest period to each value will be selected.

mtpy.utils.shapefiles_creator.create_tipper_imag_shp_from_csv(csvfile, line_length=0.03, target_epsg_code=4283)

create imagery tipper lines shape from a csv file. this function needs csv file as input. The shape is a line without arrow. Must use a GIS software such as ArcGIS to display and add an arrow at each line’s end line_length=4 how long will be the line (arrow) return: a geopandas dataframe object for further processing.

mtpy.utils.shapefiles_creator.create_tipper_real_shp_from_csv(csvfile, line_length=0.03, target_epsg_code=4283)

create tipper lines shape from a csv file. This function needs csv file as its input. The shape is a line without arrow. Must use a GIS software such as ArcGIS to display and add an arrow at each line’s end line_length=4 how long will be the line (arrow) return: a geopandas dataframe object for further processing.

mtpy.utils.shapefiles_creator.export_geopdf_to_image(geopdf, bbox, jpg_file_name, target_epsg_code=None, colorby=None, colormap=None, showfig=False)

Export a geopandas dataframe to a jpe_file, with optionally a new epsg projection. :param geopdf: a geopandas dataframe :param bbox: This param ensures that we can set a consistent display area defined by a dict with 4 keys

[MinLat, MinLon, MaxLat, MaxLon], cover all ground stations, not just this period-dependent geopdf

Parameters

• jpg_file_name (output) – path2jpeg

• target_epsg_code – 4326 etc

• showfig – If True, then display fig on screen.

Returns

mtpy.utils.shapefiles_creator.plot_phase_tensor_ellipses_and_tippers(edi_dir, out_dir, iperiod=0)

plot phase tensor ellipses and tipers into one figure. :param edi_dir: edi directory :param outfile: save figure to output file :param iperiod: the index of periods :return: saved figure file

mtpy.utils.shapefiles_creator.process_csv_folder(csv_folder, bbox_dict, target_epsg_code=4283)

process all *.csv files in a dir, ude target_epsg_code=4283 GDA94 as default. This function uses csv-files folder as its input. :param csv_folder: :return:

Create shape files for phase tensor ellipses. https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data

Created on Sun Apr 13 12:32:16 2014

@author: jrpeacock

class mtpy.utils.shapefiles.PTShapeFile(edi_list=None, proj='WGS84', esize=0.03, **kwargs)

write shape file for GIS plotting programs

key words/attributes	Description
edi_list	list of edi files, full paths
ellipse_size	size of normalized ellipse in map scale default is .01
mt_obj_list	list of mt.MT objects default is None, filled if edi_list is given
plot_period	list or value of period to convert to shape file default is None, which will write a file for every period in the edi files
ptol	tolerance to look for given periods default is .05
pt_dict	dictionary with keys of plot_period. Each dictionary key is a structured array containing the important information for the phase tensor.
projection	projection of coordinates see EPSG for all options default is WSG84 in lat and lon
save_path	path to save files to default is current working directory.

Methods	Description
_get_plot_period	get a list of all frequencies possible from input files
_get_pt_array	get phase tensors from input files and put the information into a structured array
write_shape_files	write shape files based on attributes of class

• This will project the data into UTM WSG84

Example

:: >>> edipath = r”/home/edi_files_rotated_to_geographic_north” >>> edilist = [os.path.join(edipath, edi) for edi in os.listdir(edipath) if edi.find(‘.edi’)>0] >>> pts = PTShapeFile(edilist, save_path=r”/home/gis”) >>> pts.write_shape_files()

• To project into another datum, set the projection attribute

Example

:: >>> pts = PTShapeFile(edilist, save_path=r”/home/gis”) >>> pts.projection = ‘NAD27’ >>> pts.write_shape_files()

Attributes

rotation_angle

rotation angle of Z and Tipper

Methods

write_data_pt_shape_files_modem(modem_data_fn)	write pt files from a modem data file.
write_residual_pt_shape_files_modem(...[, ...])	write residual pt shape files from ModEM output
write_resp_pt_shape_files_modem(...[, ...])	write pt files from a modem response file where ellipses are normalized by the data file.
write_shape_files([periods])	write shape file from given attributes https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html #create-a-new-shapefile-and-add-data

property rotation_angle

rotation angle of Z and Tipper

write_data_pt_shape_files_modem(modem_data_fn, rotation_angle=0.0)

write pt files from a modem data file.

write_residual_pt_shape_files_modem(modem_data_fn, modem_resp_fn, rotation_angle=0.0, normalize='1')

write residual pt shape files from ModEM output

normalize [ ‘1’ | ‘all’ ]

• ‘1’ to normalize the ellipse by itself, all ellipses are

  normalized to phimax, thus one axis is of length 1*ellipse_size

• ‘all’ to normalize each period by the largest phimax

write_resp_pt_shape_files_modem(modem_data_fn, modem_resp_fn, rotation_angle=0.0)

write pt files from a modem response file where ellipses are normalized by the data file.

write_shape_files(periods=None)

write shape file from given attributes https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html #create-a-new-shapefile-and-add-data

class mtpy.utils.shapefiles.TipperShapeFile(edi_list=None, **kwargs)

write shape file for GIS plotting programs.

currently only writes the real induction vectors.

key words/attributes	Description
arrow_direction	[ 1 | -1 ] 1 for Weise convention –> point toward conductors. default is 1 (-1 is not supported yet)
arrow_head_height	height of arrow head in map units default is .002
arrow_head_width	width of arrow head in map units default is .001
arrow_lw	width of arrow in map units default is .0005
arrow_size	size of normalized arrow length in map units default is .01
edi_list	list of edi files, full paths
mt_obj_list	list of mt.MT objects default is None, filled if edi_list is given
plot_period	list or value of period to convert to shape file default is None, which will write a file for every period in the edi files
ptol	tolerance to look for given periods default is .05
pt_dict	dictionary with keys of plot_period. Each dictionary key is a structured array containing the important information for the phase tensor.
projection	projection of coordinates see EPSG for all options default is WSG84
save_path	path to save files to default is current working directory.

Methods	Description
_get_plot_period	get a list of all possible frequencies from data
_get_tip_array	get Tipper information from data and put into a structured array for easy manipulation
write_real_shape_files	write real induction arrow shape files
write_imag_shape_files	write imaginary induction arrow shape files

Example

:: >>> edipath = r”/home/edi_files_rotated_to_geographic_north” >>> edilist = [os.path.join(edipath, edi) for edi in os.listdir(edipath) if edi.find(‘.edi’)>0] >>> tipshp = TipperShapeFile(edilist, save_path=r”/home/gis”) >>> tipshp.arrow_head_height = .005 >>> tipshp.arrow_lw = .0001 >>> tipshp.arrow_size = .05 >>> tipshp.write_shape_files()

Attributes

rotation_angle

rotation angle of Z and Tipper

Methods

write_imag_shape_files()	write shape file from given attributes
write_real_shape_files()	write shape file from given attributes
write_tip_shape_files_modem(modem_data_fn[, ...])	write tip files from a modem data file.
write_tip_shape_files_modem_residual(...)	write residual tipper files for modem

property rotation_angle

rotation angle of Z and Tipper

write_imag_shape_files()

write shape file from given attributes

write_real_shape_files()

write shape file from given attributes

write_tip_shape_files_modem(modem_data_fn, rotation_angle=0.0)

write tip files from a modem data file.

write_tip_shape_files_modem_residual(modem_data_fn, modem_resp_fn, rotation_angle)

write residual tipper files for modem

mtpy.utils.shapefiles.create_phase_tensor_shpfiles(edi_dir, save_dir, proj='WGS84', ellipse_size=1000, every_site=1, period_list=None)

generate shape file for a folder of edi files, and save the shape files a dir. :param edi_dir: :param save_dir: :param proj: defult is WGS84-UTM, with ellipse_size=1000 meters :param ellipse_size: the size of ellipse: 100-5000, try them out to suit your needs :param every_site: by default every MT station will be output, but user can sample down with 2, 3,.. :return:

mtpy.utils.shapefiles.create_tipper_shpfiles(edipath, save_dir)

Create Tipper (induction arrows real and imaginary) shape files :param edipath: :param save_dir: :return:

mtpy.utils.shapefiles.modem_to_shapefiles(mfndat, save_dir)

create shape file representaiotn for ModEM model :param mfndat: path2Modular_NLCG_110.dat :param save_dir: path2outshp :return:

mtpy.utils.shapefiles.reproject_layer(in_shape_file, out_shape_file=None, out_proj='WGS84')

reproject coordinates into a different coordinate system

GIS Tools

GIS_TOOLS

This module contains tools to help project between coordinate systems. The module will first use GDAL if installed. If GDAL is not installed then pyproj is used. A test has been made for new versions of GDAL which swap the input lat and lon when using transferPoint, so the user should not have to worry about which version they have.

Main functions are:

• project_point_ll2utm

• project_point_utm2ll

These can take in a point or an array or list of points to project.

latitude and longitude can be input as:

• ‘DD:mm:ss.ms’

• ‘DD.decimal_degrees’

• float(DD.decimal_degrees)

Created on Fri Apr 14 14:47:48 2017 Revised: 5/2020 JP

@author: jrpeacock

exception mtpy.utils.gis_tools.GISError

mtpy.utils.gis_tools.assert_elevation_value(elevation)

make sure elevation is a floating point number

mtpy.utils.gis_tools.assert_lat_value(latitude)

make sure latitude is in decimal degrees

mtpy.utils.gis_tools.assert_lon_value(longitude)

make sure longitude is in decimal degrees

mtpy.utils.gis_tools.convert_position_float2str(position)

convert position float to a string in the format of DD:MM:SS

Parameters

position (float) – decimal degrees of latitude or longitude

Return type

float

Returns

latitude or longitude in DD:MM.SS.ms

Example

:: >>> import mtpy.utils.gis_tools as gis_tools >>> gis_tools.convert_position_float2str(-118.34563) ‘-118:34:56.30’

mtpy.utils.gis_tools.convert_position_str2float(position_str)

Convert a position string in the format of DD:MM:SS to decimal degrees

Parameters

position_str (string [ 'DD:MM:SS.ms' | 'DD.degrees' ]) – degrees of latitude or longitude

Return type

float

Returns

latitude or longitude in decimal degrees

Example

>>> from mtpy.utils import gis_tools
>>> gis_tools.convert_position_str2float('-118:34:56.3')

-118.58230555555555

mtpy.utils.gis_tools.epsg_project(x, y, epsg_from, epsg_to, proj_str=None)

project some xy points using the pyproj modules

Parameters

xinteger or float

x coordinate of point

yinteger or float

y coordinate of point

epsg_fromint

epsg code of x, y points provided. To provide custom projection, set to 0 and provide proj_str

epsg_toTYPE

epsg code to project to. To provide custom projection, set to 0 and provide proj_str

proj_strstr

Proj4 string to provide to pyproj if using custom projection. This proj string will be applied if epsg_from or epsg_to == 0. The default is None.

Returns

xp, yp

x and y coordinates of projected point.

mtpy.utils.gis_tools.get_epsg(latitude, longitude)

get epsg code for the utm projection (WGS84 datum) of a given latitude and longitude pair

Parameters

• latitude ([ string | float ]) – latitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

• longitude ([ string | float ]) – longitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

Returns

EPSG number

Return type

int

Example

>>> gis_tools.get_epsg(-34.299442, '149:12:03.71')

32755

mtpy.utils.gis_tools.get_utm_zone(latitude, longitude)

Get utm zone from a given latitude and longitude

Parameters

• latitude ([ string | float ]) – latitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

• longitude ([ string | float ]) – longitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

Returns

zone number

Return type

int

Returns

is northern

Return type

[ True | False ]

Returns

UTM zone

Return type

string

Example

>>> lat, lon = ('-34:17:57.99', 149.2010301)
>>> zone_number, is_northing, utm_zone = gis_tools.get_utm_zone(lat, lon)
>>> print(zone_number, is_northing, utm_zone)

(55, False, ‘55H’)

mtpy.utils.gis_tools.project_point_ll2utm(lat, lon, datum='WGS84', utm_zone=None, epsg=None)

Project a point that is in latitude and longitude to the specified UTM coordinate system.

Parameters

• latitude ([ string | float ]) – latitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

• longitude ([ string | float ]) – longitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

• datum (string) – well known datum

• utm_zone ([ string | int ]) – utm_zone {0-9}{0-9}{C-X} or {+, -}{0-9}{0-9}

• epsg ([ int | string ]) – EPSG number defining projection (see http://spatialreference.org/ref/ for moreinfo) Overrides utm_zone if both are provided

Returns

project point(s)

Return type

tuple if a single point, np.recarray if multiple points * tuple is (easting, northing,utm_zone) * recarray has attributes (easting, northing, utm_zone, elevation)

Single Point

>>> gis_tools.project_point_ll2utm('-34:17:57.99', '149.2010301')

(702562.6911014864, 6202448.5654573515, ‘55H’)

Multiple Points

>>> lat = np.arange(20, 40, 5)
>>> lon = np.arange(-110, -90, 5)
>>> gis_tools.project_point_ll2utm(lat, lon, datum='NAD27')

rec.array([( -23546.69921068, 2219176.82320353, 0., ‘13R’),

( 500000. , 2764789.91224626, 0., ‘13R’), ( 982556.42985037, 3329149.98905941, 0., ‘13R’), (1414124.6019547 , 3918877.48599922, 0., ‘13R’)],

dtype=[(‘easting’, ‘<f8’), (‘northing’, ‘<f8’),

(‘elev’, ‘<f8’), (‘utm_zone’, ‘<U3’)])

mtpy.utils.gis_tools.project_point_utm2ll(easting, northing, utm_zone, datum='WGS84', epsg=None)

Project a point that is in UTM to the specified geographic coordinate system.

Parameters

• easting (float) – easting in meters

• northing (float) – northing in meters

• datum (string) – well known datum

• utm_zone ([ string | int ]) – utm_zone {0-9}{0-9}{C-X} or {+, -}{0-9}{0-9}

• epsg ([ int | string ]) – EPSG number defining projection (see http://spatialreference.org/ref/ for moreinfo) Overrides utm_zone if both are provided

Returns

project point(s)

Return type

tuple if a single point, np.recarray if multiple points * tuple is (easting, northing,utm_zone) * recarray has attributes (easting, northing, utm_zone, elevation)

Single Point

>>> gis_tools.project_point_utm2ll(670804.18810336,

… 4429474.30215206, … datum=’WGS84’, … utm_zone=’11T’, … epsg=26711) (40.000087, -114.999128)

Multiple Points

>>> gis_tools.project_point_utm2ll([670804.18810336, 680200],

… [4429474.30215206, 4330200], … datum=’WGS84’, utm_zone=’11T’, … epsg=26711) rec.array([(40.000087, -114.999128), (39.104208, -114.916058)],

dtype=[(‘latitude’, ‘<f8’), (‘longitude’, ‘<f8’)])

mtpy.utils.gis_tools.split_utm_zone(utm_zone)

Split utme zone into zone number and is northing

Parameters

utm_zone ([ string | int ]) – utm zone string as {0-9}{0-9}{C-X} or {+,-}{0-9}{0-9}

Returns

utm zone number

Return type

int

Returns

is_northern

Return type

boolean [ True | False ]

Example

>>> gis_tools.split_utm_zone('11S')

11, True

mtpy.utils.gis_tools.utm_letter_designator(latitude)

Get the UTM zone letter designation for a given latitude

Parameters

latitude ([ string | float ]) – latitude in [ ‘DD:mm:ss.ms’ | ‘DD.decimal’ | float ]

Returns

UTM zone letter designation

Return type

string

Example

>>> gis_utils.utm_letter_designator('-34:17:57.99')

H

mtpy.utils.gis_tools.utm_wgs84_conv(lat, lon)

Bidirectional UTM-WGS84 converter https://github.com/Turbo87/utm/blob/master/utm/conversion.py :param lat: :param lon: :return: tuple(e, n, zone, lett)

mtpy.utils.gis_tools.utm_zone_to_epsg(zone_number, is_northern)

get epsg code (WGS84 datum) for a given utm zone

Parameters

• zone_number (int) – UTM zone number

• is_northing ([ True | False ]) – Boolean of UTM is in northern hemisphere

Returns

EPSG number

Return type

int

Example

>>> gis_tools.utm_zone_to_epsg(55, False)

32755

mtpy.utils.gis_tools.validate_epsg(epsg)

Make sure epsg is an integer

Parameters

epsg ([ int | str ]) – EPSG number

Returns

EPSG number

Return type

int

mtpy.utils.gis_tools.validate_input_values(values, location_type=None)

make sure the input values for lat, lon, easting, northing will be an numpy array with a float data type

can input a string as a comma separated list

Parameters

values ([ float | string | list | numpy.ndarray ]) – values to project, can be given as: * float * string of a single value or a comma separate string ‘34.2, 34.5’ * list of floats or string * numpy.ndarray

Returns

array of floats

Return type

numpy.ndarray(dtype=float)

mtpy.utils.gis_tools.validate_utm_zone(utm_zone)

Make sure utm zone is a valid string

Parameters

utm_zone ([ int | string ]) – UTM zone as {0-9}{0-9}{C-X} or {+, -}{0-9}{0-9}

Returns

valid UTM zone

Return type

[ int | string ]

Other Tools

Created on Wed Oct 25 09:35:31 2017

@author: Alison Kirkby

functions to assist with mesh generation

mtpy.utils.mesh_tools.get_nearest_index(array, value)

Return the index of the nearest value to the provided value in an array:

inputs:

array = array or list of values value = target value

mtpy.utils.mesh_tools.get_padding_cells(cell_width, max_distance, num_cells, stretch)

get padding cells, which are exponentially increasing to a given distance. Make sure that each cell is larger than the one previously.

Returns

**padding**np.ndarray

array of padding cells for one side

mtpy.utils.mesh_tools.get_padding_cells2(cell_width, core_max, max_distance, num_cells)

get padding cells, which are exponentially increasing to a given distance. Make sure that each cell is larger than the one previously.

mtpy.utils.mesh_tools.get_padding_from_stretch(cell_width, pad_stretch, num_cells)

get padding cells using pad stretch factor

mtpy.utils.mesh_tools.get_rounding(cell_width)

Get the rounding number given the cell width. Will be one significant number less than the cell width. This reduces weird looking meshes.

Parameters

cell_width (float) – Width of mesh cell

Returns

digit to round to

Return type

int

>>> from mtpy.utils.mesh_tools import get_rounding
>>> get_rounding(9)
0
>>> get_rounding(90)
-1
>>> get_rounding(900)
-2
>>> get_rounding(9000)
-3

mtpy.utils.mesh_tools.get_station_buffer(grid_east, grid_north, station_east, station_north, buf=10000.0)

get cells within a specified distance (buf) of the stations returns a 2D boolean (True/False) array

mtpy.utils.mesh_tools.grid_centre(grid_edges)

calculate the grid centres from an array that defines grid edges :param grid_edges: array containing grid edges :returns: grid_centre: centre points of grid

mtpy.utils.mesh_tools.interpolate_elevation_to_grid(grid_east, grid_north, epsg=None, utm_zone=None, surfacefile=None, surface=None, method='linear', fast=True, buffer=1)

# Note: this documentation is outdated and seems to be copied from # model.interpolate_elevation2. It needs to be updated. This # funciton does not update a dictionary but returns an array of # elevation data.

project a surface to the model grid and add resulting elevation data to a dictionary called surface_dict. Assumes the surface is in lat/long coordinates (wgs84) The ‘fast’ method extracts a subset of the elevation data that falls within the mesh-bounds and interpolates them onto mesh nodes. This approach significantly speeds up (~ x5) the interpolation procedure.

returns nothing returned, but surface data are added to surface_dict under the key given by surfacename.

inputs choose to provide either surface_file (path to file) or surface (tuple). If both are provided then surface tuple takes priority.

surface elevations are positive up, and relative to sea level. surface file format is:

ncols 3601 nrows 3601 xllcorner -119.00013888889 (longitude of lower left) yllcorner 36.999861111111 (latitude of lower left) cellsize 0.00027777777777778 NODATA_value -9999 elevation data W –> E N | V S

Alternatively, provide a tuple with: (lon,lat,elevation) where elevation is a 2D array (shape (ny,nx)) containing elevation points (order S -> N, W -> E) and lon, lat are either 1D arrays containing list of longitudes and latitudes (in the case of a regular grid) or 2D arrays with same shape as elevation array containing longitude and latitude of each point.

other inputs: surfacename = name of surface for putting into dictionary surface_epsg = epsg number of input surface, default is 4326 for lat/lon(wgs84) method = interpolation method. Default is ‘nearest’, if model grid is dense compared to surface points then choose ‘linear’ or ‘cubic’

mtpy.utils.mesh_tools.make_log_increasing_array(z1_layer, target_depth, n_layers, increment_factor=0.9)

create depth array with log increasing cells, down to target depth, inputs are z1_layer thickness, target depth, number of layers (n_layers)

mtpy.utils.mesh_tools.rotate_mesh(grid_east, grid_north, origin, rotation_angle, return_centre=False)

rotate a mesh defined by grid_east and grid_north.

Parameters

• grid_east – 1d array defining the edges of the mesh in the east-west direction

• grid_north – 1d array defining the edges of the mesh in the north-south direction

• origin – real-world position of the (0,0) point in grid_east, grid_north

• rotation_angle – angle in degrees to rotate the grid by

• return_centre – True/False option to return points on centre of grid instead of grid edges

Returns

grid_east, grid_north - 2d arrays describing the east and north coordinates

A more Pythonic way of logging: Define a class MtPyLog to wrap the python logging module; Use a (optional) configuration file (yaml, ini, json) to configure the logging, It will return a logger object with the user-provided config setting. see also: http://www.cdotson.com/2015/11/python-logging-best-practices/