2015-11-11 9 views
7
का उपयोग कर

enter image description hereकी साजिश रचने भूकंपीय लचीलेपन की matplotlib

मैं matplotlib का उपयोग कर की साजिश रचने के ऊपर शैली पुन: बनाने की कोशिश कर रहा हूँ निशान बनता है।

कच्चा डेटा 2 डी numpy सरणी में संग्रहीत किया जाता है, जहां तेज अक्ष समय होता है।

लाइनों को प्लॉट करना आसान बिट है। मैं छायांकित क्षेत्रों में कुशलतापूर्वक प्राप्त करने की कोशिश कर रहा हूं।

मेरे वर्तमान प्रयास तरह दिखता है:

import numpy as np 
from matplotlib import collections 
import matplotlib.pyplot as pylab 

#make some oscillating data 
panel = np.meshgrid(np.arange(1501), np.arange(284))[0] 
panel = np.sin(panel) 

#generate coordinate vectors. 
panel[:,-1] = np.nan #lazy prevents polygon wrapping 
x = panel.ravel() 
y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel() 

#find indexes of each zero crossing 
zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1 

#calculate scalars used to shift "traces" to plotting corrdinates 
trace_centers = np.linspace(1,284, panel.shape[-2]).reshape(-1,1) 
gain = 0.5 #scale traces 

#shift traces to plotting coordinates 
x = ((panel*gain)+trace_centers).ravel() 

#split coordinate vectors at each zero crossing 
xpoly = np.split(x, zero_crossings) 
ypoly = np.split(y, zero_crossings) 

#we only want the polygons which outline positive values 
if x[0] > 0: 
    steps = range(0, len(xpoly),2) 
else: 
    steps = range(1, len(xpoly),2) 

#turn vectors of polygon coordinates into lists of coordinate pairs 
polygons = [zip(xpoly[i], ypoly[i]) for i in steps if len(xpoly[i]) > 2] 

#this is so we can plot the lines as well 
xlines = np.split(x, 284) 
ylines = np.split(y, 284) 
lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))] 

#and plot 
fig = pylab.figure() 
ax = fig.add_subplot(111) 
col = collections.PolyCollection(polygons) 
col.set_color('k') 
ax.add_collection(col, autolim=True) 
col1 = collections.LineCollection(lines) 
col1.set_color('k') 
ax.add_collection(col1, autolim=True) 
ax.autoscale_view() 
pylab.xlim([0,284]) 
pylab.ylim([0,1500]) 
ax.set_ylim(ax.get_ylim()[::-1]) 
pylab.tight_layout() 
pylab.show() 

और परिणाम enter image description here

है दो मुद्दों के होते हैं:

  1. यह पूरी तरह से भर नहीं है क्योंकि मैं पर बंटवारे हूँ शून्य क्रॉसिंग के निकटतम सरणी इंडेक्स, सटीक शून्य क्रॉसिंग नहीं। मुझे लगता है कि प्रत्येक शून्य क्रॉसिंग की गणना एक बड़ी कम्प्यूटेशनल हिट होगी।

  2. प्रदर्शन। यह बुरा नहीं है, समस्या के आकार को देखते हुए - मेरे लैपटॉप पर रेंडर करने के लिए लगभग एक सेकंड, लेकिन मैं इसे 100ms - 200ms तक ले जाना चाहता हूं।

उपयोग के मामले में मैं अजगर/scipy/matplotlib के साथ अजगर तक सीमित हूं। कोई सुझाव?

फॉलोअप:

रैखिक शून्य क्रॉसिंग interpolating बाहर कर देता है बहुत कम कम्प्यूटेशनल लोड के साथ किया जा सकता है। डेटा में इंटरपोलेटेड मान डालने, नैनों को नकारात्मक मान सेट करके और pyplot.fill पर एकल कॉल का उपयोग करके, 500,000 विषम नमूने लगभग 300ms में प्लॉट किए जा सकते हैं।

संदर्भ के लिए, उसी डेटा पर नीचे टॉम की विधि लगभग 8 सेकंड लग गई।

निम्न कोड एक डिस्प्ले के साथ एक numpy recarray का इनपुट मानता है जो एक भूकंपीय यूनिक्स शीर्षलेख/ट्रेस परिभाषा की नकल करता है।

def wiggle(frame, scale=1.0): 
     fig = pylab.figure() 
     ax = fig.add_subplot(111)   
     ns = frame['ns'][0] 
     nt = frame.size 
     scalar = scale*frame.size/(frame.size*0.2) #scales the trace amplitudes relative to the number of traces 
     frame['trace'][:,-1] = np.nan #set the very last value to nan. this is a lazy way to prevent wrapping 
     vals = frame['trace'].ravel() #flat view of the 2d array. 
     vect = np.arange(vals.size).astype(np.float) #flat index array, for correctly locating zero crossings in the flat view 
     crossing = np.where(np.diff(np.signbit(vals)))[0] #index before zero crossing 
     #use linear interpolation to find the zero crossing, i.e. y = mx + c. 
     x1= vals[crossing] 
     x2 = vals[crossing+1] 
     y1 = vect[crossing] 
     y2 = vect[crossing+1] 
     m = (y2 - y1)/(x2-x1) 
     c = y1 - m*x1  
     #tack these values onto the end of the existing data 
     x = np.hstack([vals, np.zeros_like(c)]) 
     y = np.hstack([vect, c]) 
     #resort the data 
     order = np.argsort(y) 
     #shift from amplitudes to plotting coordinates 
     x_shift, y = y[order].__divmod__(ns) 
     ax.plot(x[order] *scalar + x_shift + 1, y, 'k') 
     x[x<0] = np.nan 
     x = x[order] *scalar + x_shift + 1 
     ax.fill(x,y, 'k', aa=True) 
     ax.set_xlim([0,nt]) 
     ax.set_ylim([ns,0]) 
     pylab.tight_layout() 
     pylab.show() 

enter image description here

पूर्ण कोड https://github.com/stuliveshere/PySeis

उत्तर

5

पर प्रकाशित किया जाता है आप fill_betweenx साथ आसानी से कर सकते हैं। दस्तावेज़ों से:

दो क्षैतिज घटता के बीच बहुभुज बनाओ।

कॉल हस्ताक्षर:

fill_betweenx (y, x1, x2 = 0, जहां कोई नहीं =, ** kwargs) एक PolyCollection x1 और x2 के बीच क्षेत्रों भरने बनाएं जहां जहां == सच

यहां महत्वपूर्ण हिस्सा where तर्क है।

तो, आप x2 = offset करना चाहते हैं, और उसके बाद where = x>offset

उदाहरण के लिए है:

import numpy as np 
import matplotlib.pyplot as plt 

fig,ax = plt.subplots() 

# Some example data 
y = np.linspace(700.,900.,401) 
offset = 94. 
x = offset+10*(np.sin(y/2.)* 
     1/(10. * np.sqrt(2 * np.pi)) * 
     np.exp(- (y - 800)**2/(2 * 10.**2)) 
     ) # This function just gives a wave that looks something like a seismic arrival 

ax.plot(x,y,'k-') 
ax.fill_betweenx(y,offset,x,where=(x>offset),color='k') 

ax.set_xlim(93,95) 

plt.show() 

enter image description here

आप अपने ऑफसेट से प्रत्येक के लिए fill_betweenx करने की ज़रूरत है। उदाहरण के लिए:

import numpy as np 
import matplotlib.pyplot as plt 

fig,ax = plt.subplots() 

# Some example data 
y = np.linspace(700.,900.,401) 
offsets = [94., 95., 96., 97.] 
times = [800., 790., 780., 770.] 

for offset, time in zip(offsets,times): 
    x = offset+10*(np.sin(y/2.)* 
     1/(10. * np.sqrt(2 * np.pi)) * 
     np.exp(- (y - time)**2/(2 * 10.**2)) 
     ) 

    ax.plot(x,y,'k-') 
    ax.fill_betweenx(y,offset,x,where=(x>offset),color='k') 

ax.set_xlim(93,98) 

plt.show() 

enter image description here

+0

भूखंडों खूबसूरती से, लेकिन रूपरेखा इसके चारों ओर 5 बार मेरी विधि की तुलना में धीमी है पता चलता है, मैं क्योंकि आप प्रत्येक का पता लगाने के ऊपर पुनरावृति करने के लिए है यह मानते हुए कि कर रहा हूँ, ताकि आप की साजिश रचने कर रहे हैं एक बड़े से अधिक सैकड़ों छोटे संग्रह। मैं आज रात थोड़ा गहरा प्रोफाइलिंग में खोदना होगा। – scrooge

1

यह काफी अगर तुम SEGY प्रारूप और/या txt प्रारूप में अपने भूकंपीय निशान है करने के लिए आसान है (आप .txt प्रारूप में उन्हें अंत में की आवश्यकता होगी)। सबसे अच्छी विधि खोजने में काफी समय लगा। पाइथन और प्रोग्रामिंग के लिए भी नया नया, तो कृपया नम्र रहें।

सेजी फ़ाइल को एक .txt फ़ाइल में कनवर्ट करने के लिए मैंने सेईसी (http://dmng.ru/en/freeware.html; रूसी साइट पर ध्यान न दें, यह एक वैध प्रोग्राम है)। लोड करने और प्रदर्शित करने के लिए आपको numpy और matplotlib की आवश्यकता है।

निम्नलिखित कोड भूकंपीय निशान लोड करेंगे, उन्हें स्थानांतरित करेंगे, और उन्हें साजिश करेंगे। जाहिर है आपको अपनी फाइल लोड करने, लंबवत और क्षैतिज श्रेणियों को बदलने की जरूरत है, और vmin और vmax के साथ थोड़ा सा खेलना होगा। यह एक भूरे रंग के रंग का भी उपयोग करता है। कोड इस तरह एक छवि का उत्पादन करेगा: http://goo.gl/0meLyz

import numpy as np 
import matplotlib.pyplot as plt 

traces = np.loadtxt('yourtracestxtfile.txt') 
traces = np.transpose(traces) 

seismicplot = plt.imshow(traces[3500:4500,500:900], cmap = 'Greys',vmin = 0,vmax = 1,aspect = 'auto') #Tip: traces[vertical range,horizontal range] 
संबंधित मुद्दे