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

class mtpy.utils.shapefiles_creator.ShapeFilesCreator(edifile_list, outdir, orig_crs={'init': 'epsg:4283'})

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

create_measurement_csv(self, 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(self[, outshpfile])	create station location geopandas dataframe, and output to shape file
create_phase_tensor_csv(self, dest_dir[, …])	create phase tensor ellipse and tipper properties.
create_phase_tensor_csv_with_image(\*args, …)	Using PlotPhaseTensorMaps class to generate csv file of phase tensor attributes, etc.
create_phase_tensor_shp(self, period[, …])	create phase tensor ellipses shape file correspond to a MT period :return: (geopdf_obj, path_to_shapefile)
create_tipper_imag_shp(self, period[, …])	create imagery tipper lines shapefile from a csv file The shapefile consists of lines without arrow.
create_tipper_real_shp(self, period[, …])	create real tipper lines shapefile from a csv file The shapefile consists of lines without arrow.
display_on_basemap(self)	display MT stations which are in stored in geopandas dataframe in a base map.
display_on_image(self)	display/overlay the MT properties on a background geo-referenced map image
export_edi_files(self, dest_dir[, …])	export edi files.
get_bounding_box(self[, epsgcode])	compute bounding box
get_min_max_distance(self)	get the min and max distance between all possible pairs of stations.
get_period_occurance(self, 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(self[, 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(self, period[, …])	For a given MT period (s) value, compute the phase tensor and tippers etc.
get_station_utmzones_stats(self)	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(self)	get the min max statistics of the distances between stations.
plot_stations(self[, savefile, showfig])	Visualise the geopandas df of MT stations
select_periods(self[, show, period_list, …])	Use edi_collection to analyse the whole set of EDI files
show_obj(self[, dest_dir])	test call object’s methods and show it’s properties

create_phase_tensor_shp(self, 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(self, 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(self, 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_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, outfile=None, 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(self, …[, …])	write pt files from a modem data file.
write_residual_pt_shape_files_modem(self, …)	write residual pt shape files from ModEM output
write_resp_pt_shape_files_modem(self, …[, …])	write pt files from a modem response file where ellipses are normalized by the data file.
write_shape_files(self[, periods])	write shape file from given attributes https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html #create-a-new-shapefile-and-add-data

rotation_angle

rotation angle of Z and Tipper

write_data_pt_shape_files_modem(self, modem_data_fn, rotation_angle=0.0)

write pt files from a modem data file.

write_residual_pt_shape_files_modem(self, 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(self, 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(self, 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(self)	write shape file from given attributes
write_real_shape_files(self)	write shape file from given attributes
write_tip_shape_files_modem(self, modem_data_fn)	write tip files from a modem data file.
write_tip_shape_files_modem_residual(self, …)	write residual tipper files for modem

rotation_angle

rotation angle of Z and Tipper

write_imag_shape_files(self)

write shape file from given attributes

write_real_shape_files(self)

write shape file from given attributes

write_tip_shape_files_modem(self, modem_data_fn, rotation_angle=0.0)

write tip files from a modem data file.

write_tip_shape_files_modem_residual(self, 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

Created on Fri Apr 14 14:47:48 2017

@author: jrpeacock

exception mtpy.utils.gis_tools.GIS_ERROR

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

Returns:

**position_str** : string

latitude or longitude in format of DD:MM:SS.ms

mtpy.utils.gis_tools.convert_position_str2float(position_str)

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

Returns:

**position** : float

latitude or longitude in decimal degrees

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

project some xy points using the pyproj modules

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

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

mtpy.utils.gis_tools.get_utm_string_from_sr(*args, **kwargs)

return utm zone string from spatial reference instance

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

Get utm zone from a given latitude and longitude

mtpy.utils.gis_tools.ll_to_utm(*args, **kwargs)

converts lat/long to UTM coords. Equations from USGS Bulletin 1532 East Longitudes are positive, West longitudes are negative. North latitudes are positive, South latitudes are negative Lat and Long are in decimal degrees Written by Chuck Gantz- chuck.gantz@globalstar.com

Outputs:

UTMzone, easting, northing

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

Project a point that is in Lat, Lon (will be converted to decimal degrees) into UTM coordinates.

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

Project a point that is in Lat, Lon (will be converted to decimal degrees) into UTM coordinates.

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

Project a list of points that is in Lat, Lon (will be converted to decimal degrees) into UTM coordinates.

mtpy.utils.gis_tools.utm_to_ll(*args, **kwargs)

converts UTM coords to lat/long. Equations from USGS Bulletin 1532 East Longitudes are positive, West longitudes are negative. North latitudes are positive, South latitudes are negative Lat and Long are in decimal degrees. Written by Chuck Gantz- chuck.gantz@globalstar.com Converted to Python by Russ Nelson <nelson@crynwr.com>

Outputs:

Lat,Lon

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

Other Tools

class mtpy.utils.decorator.deprecated(reason)

Description:

used to mark functions, methods and classes deprecated, and prints warning message when it called decorators based on https://stackoverflow.com/a/40301488

Usage:

todo: write usage

Author: YingzhiGou Date: 20/06/2017

Methods

__call__

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_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)

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/