Fei Zhang
2018-10-23
In [1]:
import os
import sys
import glob
import fnmatch
from mtpy.utils.shapefiles_creator import ShapeFilesCreator
In [2]:
def recursive_glob(dirname, ext='*.edi'):
"""
Under the dirname recursively find all files with extension ext.
Return a list of the full-path to the types of files of interest.
This function is useful to handle a nested directories of EDI files.
:param dirname: a single dir OR a list of dirs.
:param ext: eg, ".xml"
:return: a list of path2files
"""
if isinstance(dirname, (list,)): # the input argument is a list of directory
filelist = []
for adir in dirname:
filelist.extend(recursive_glob(adir))
return filelist
else: # input variable is a single dir
matches = []
for root, dirnames, filenames in os.walk(dirname):
for filename in fnmatch.filter(filenames, ext):
matches.append(os.path.join(root, filename))
return matches
In [3]:
EDI_DIR=r'C:\mtpywin\mtpy\examples\data\edi_files_2'
EDI_DIR=r'C:\mtpywin\mtpy\examples\data\edi2'
# EDI_DIR=r'C:\mtpywin\mtpy\examples\data\edi_files' # no tipper
# EDI_DIR = 'E:/Data/MT_Datasets/3D_MT_data_edited_fromDuanJM'
# EDI_DIR = r"E:\Data\MT_Datasets\GA_UA_edited_10s-10000s"
# EDI_DIR = r"E:\Data\MT_Datasets\Isa_EDI_edited_10Hz_1000s"
# EDI_DIR = r"E:\Data\MT_Datasets\728889\EDI_files" # narrow area, not shown
# EDI_DIR = r"E:\Data\MT_Datasets\75098\EDI_files" # cross UTM zones 51 and 52, No tippers?
# EDI_DIR =r"E:\Data\MT_Datasets\75099_Youanmi\EDI_Files_edited"
edifiles=recursive_glob(EDI_DIR)
print("Number of EDI files found = %s"%len(edifiles))
In [4]:
shpobj = ShapeFilesCreator(edifiles, "c:/tmp20181023")
In [5]:
allper=shpobj.all_unique_periods
i=10 # which period index?
In [18]:
# set your projection, this can be a UTM zone projection such as 32754 = S54
my_proj_epsg = 32754 #None # 3112 # projection GDA94/GALCC =3112
gpd_phtensor = shpobj.create_phase_tensor_shp(allper[i],target_epsg_code= my_proj_epsg)[0]
In [19]:
type(gpd_phtensor)
Out[19]:
In [20]:
gpd_phtensor.plot(linewidth=2, facecolor='grey', edgecolor='black')
Out[20]:
In [21]:
gpd_retip = shpobj.create_tipper_real_shp(allper[i],target_epsg_code= my_proj_epsg)[0]
gpd_imtip = shpobj.create_tipper_imag_shp(allper[i],target_epsg_code= my_proj_epsg)[0]
In [22]:
gpd_retip.plot()
gpd_imtip.plot()
Out[22]:
In [23]:
gpd_retip.info()
In [24]:
gpd_retip.head()
Out[24]:
In [25]:
import matplotlib.pyplot as plt
import geopandas as gpd
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
In [26]:
# composing two layers in a map
f, ax = plt.subplots(1, figsize=(20, 12))
# ax.set_xlim([140.5,141])
# ax.set_ylim([-21,-20])
# Add layer of polygons on the axis
gpd_phtensor.plot(ax=ax, linewidth=2, facecolor='grey', edgecolor='black')
gpd_retip.plot(ax=ax, color='red', linewidth=4)
gpd_imtip.plot(ax=ax, color='blue', linewidth=4)
# world.plot(ax=ax, alpha=0.5) # background map
plt.savefig('phasetensor_tippers.png')
plt.axis('equal')
# Display
#plt.show()
Out[26]:
In [27]:
print(gpd_imtip['tip_mag_im'].min(),gpd_imtip['tip_mag_im'].max())
In [28]:
print(gpd_retip['tip_mag_re'].min(), gpd_retip['tip_mag_re'].max())
In [29]:
shpobj.stations_distances
Out[29]:
In [ ]:
In [ ]: