diff --git a/MyProgressMonitorNew.m b/MyProgressMonitorNew.m
new file mode 100644
index 00000000..207d3c25
--- /dev/null
+++ b/MyProgressMonitorNew.m
@@ -0,0 +1,102 @@
+ classdef MyProgressMonitorNew < matlab.net.http.ProgressMonitor
+ properties
+ ProgHandle
+ Direction matlab.net.http.MessageType
+ Value uint64
+ NewDir matlab.net.http.MessageType = matlab.net.http.MessageType.Request
+ end
+
+ methods
+ function obj = MyProgressMonitorNew
+ obj.Interval = .01;
+ end
+
+ function done(obj)
+ obj.closeit();
+ end
+
+ function delete(obj)
+ obj.closeit();
+ end
+
+ function set.Direction(obj, dir)
+ obj.Direction = dir;
+ obj.changeDir();
+ end
+
+ function set.Value(obj, value)
+ obj.Value = value;
+ obj.update();
+ end
+ end
+
+ methods (Access = private)
+ function update(obj,~)
+ % called when Value is set
+ import matlab.net.http.*
+ if ~isempty(obj.Value)
+ if isempty(obj.Max)
+ % no maximum means we don't know length, so message changes on
+ % every call
+ value = 0;
+ if obj.Direction == MessageType.Request
+ msg = sprintf('Sent %d bytes...', obj.Value);
+ else
+ msg = sprintf('Received %d bytes...', obj.Value);
+ end
+ else
+ % maximum known; update proportional value
+ value = double(obj.Value)/double(obj.Max);
+ if obj.NewDir == MessageType.Request
+ % message changes only on change of direction
+ if obj.Direction == MessageType.Request
+ msg = 'Sending...';
+ else
+ msg = 'Receiving...';
+ end
+ end
+ end
+ if isempty(obj.ProgHandle)
+ % if we don't have a progress bar, display it for first time
+ obj.ProgHandle = ...
+ waitbar(value, msg, 'CreateCancelBtn', @(~,~)cancelAndClose(obj));
+
+ obj.NewDir = MessageType.Response;
+ elseif obj.NewDir == MessageType.Request || isempty(obj.Max)
+ % on change of direction or if no maximum known, change message
+ waitbar(value, obj.ProgHandle, msg);
+ obj.NewDir = MessageType.Response;
+ else
+ % no direction change else just update proportional value
+ waitbar(value, obj.ProgHandle);
+ end
+ end
+
+ function cancelAndClose(obj)
+ % Call the required CancelFcn and then close our progress bar. This is
+ % called when user clicks cancel or closes the window.
+ obj.CancelFcn();
+ obj.closeit();
+ end
+ end
+
+ function changeDir(obj,~)
+ % Called when Direction is set or changed. Leave the progress bar displayed.
+ obj.NewDir = matlab.net.http.MessageType.Request;
+ end
+ end
+
+ methods (Access=private)
+ function closeit(obj)
+ % Close the progress bar by deleting the handle so CloseRequestFcn isn't
+ % called, because waitbar calls our cancelAndClose(), which would cause
+ % recursion.
+ if ~isempty(obj.ProgHandle)
+ delete(obj.ProgHandle);
+ obj.ProgHandle = [];
+ end
+ end
+ end
+ end
+
+
\ No newline at end of file
diff --git a/importfilters/import_fdsn_event.m b/importfilters/import_fdsn_event.m
index f91c9ebf..32a4642e 100644
--- a/importfilters/import_fdsn_event.m
+++ b/importfilters/import_fdsn_event.m
@@ -127,125 +127,58 @@
disp(['sending request to:' baseurl 'query with options'])
disp(varargin)
- hasParallel=license('test','Distrib_Computing_Toolbox') && ZG.ParallelProcessingOpts.Enable;
- try
- if hasParallel
- fn = tempname;
- p=gcp();
- f=parfeval(p,@websave,1,fn, [baseurl 'query'], varargin{:},'format','text',options);
- warning('off','MATLAB:table:ModifiedAndSavedVarnames');
- nQuakes = 0;
- lastsize = 0;
- infmt = string(missing);
- earliest = datetime(missing);
- latest=datetime(missing);
- myfig = figure('Name','Downloading events');
- ax1=subplot(2,1,1);
- scDate=scatter(ax1,NaT,nan,'.');title(ax1,'Date of downloaded events');
- ax2=subplot(2,1,2);
- scLoc=scatter(ax2,nan,nan,'.');title(ax2,'Location of downloaded events');
- while f.State=="running"
- startDownloadTime = f.StartDateTime;
- pause(1);
- d=dir(fn);
- if isempty(d) || d.bytes == lastsize
- fprintf('.');
- continue
- end
- lastsize=d.bytes;
- rightnow = datetime;
- rightnow.TimeZone = f.StartDateTime.TimeZone;
- elapsedTime = rightnow - f.StartDateTime;
- bytes_per_second = lastsize / seconds(elapsedTime);
- try
- tb=readtable(fn);
- catch ME
- continue
- end
- if ~isempty(tb) && any([ismissing(tb.Time(end)) ismissing(tb.Magnitude(end))])
- tb(end,:)=[];
- end
- if isempty(tb)
- continue
- end
- events_per_minute = height(tb) / minutes(elapsedTime);
- if ismissing(infmt)
- infmt = "uuuu-MM-dd'T'HH:mm:ss";
- infmt{1}([5 8])=tb.Time{1}(5);
- infmt{1}(12)=tb.Time{1}(11);
- end
- tb.Time = datetime(extractBefore(tb.Time,20),'InputFormat',infmt);
- if isvalid(scDate)
- set(scDate,'XData',tb.Time,'YData',tb.Magnitude);
- end
- if isvalid(scLoc)
- set(scLoc,'XData',tb.Longitude,'YData',tb.Latitude);
- end
- drawnow nocallbacks
- nQuakes=height(tb);
- minTime = min(tb.Time);
- maxTime=max(tb.Time);
- minMag=min(tb.Magnitude);
- maxMag=max(tb.Magnitude);
- if minTime ~= earliest
- earliest=minTime;
- earlyStr = ""+string(earliest)+"";
- else
- earlyStr=string(earliest);
- end
- if maxTime ~= latest
- latest=maxTime;
- lateStr = ""+string(earliest)+"";
- else
- lateStr=string(latest);
- end
- timespan_per_minute = (latest - earliest) / minutes(elapsedTime);
- if timespan_per_minute > years(1)
- timespan_per_minute.Format='y';
- elseif timespan_per_minute > days(5)
- timespan_per_minute.Format='d';
- end
-
- msg.dbfprintf("[%d events found... still downloading ]\n"+...
- " Download Statistics:\n"+...
- " Elapsed Time : %s\n"+...
- " Bytes / sec : %d\n"+...
- " Events / min : %d\n"+...
- " Timespan / min : %s\n\n" +...
- " Catalog Statistics:\n"+...
- " start time : %s\n"+...
- " end time : %s\n"+...
- " magnitude range : [%g to %g]\n",...
- nQuakes, elapsedTime,...
- round(bytes_per_second), round(events_per_minute), string(timespan_per_minute),...
- earlyStr, lateStr, min(tb.Magnitude), max(tb.Magnitude));
- fprintf('\nwaiting for next update');
+ % hasParallel=license('test','Distrib_Computing_Toolbox') && ZG.ParallelProcessingOpts.Enable;
+ myuri=[baseurl, 'query?', sprintf('%s=%s&',string(varargin)), 'format=text'];
+ resp = get_low_level_fdsn_query(myuri);
+ switch resp.StatusCode
+ case "OK"
+ ok=true;
+ case "NoContent"
+ warning("No Data was found");
+ uOutput=[];
+ ok=false;
+ case "BadRequest"
+ disp(resp.Body.Data)
+
+ % as of 2018-12-14, USGS returns this result when limit is exceeded. these depend on the error message wording
+ maxSearchLimit = double(extractBefore(extractAfter(resp.Body.Data,'exceeds search limit of '),'.'));
+ nFound = double(extractBefore(extractAfter(resp.Body.Data,'exceeds search limit of '),'.'));
+ if ~ismissing(maxSearchLimit)
+ warning("maximum number of events [%d] exceeded. atttempting to limit results", maxSearchLimit);
+ % try again, in chunks of maxSearchLimit
+ disp('* trying again while limiting results')
+ [resp, ok] = get_in_chunks(myuri, maxSearchLimit, nFound);
+ else
+ uOutput=[];
+ ok=false;
end
- msg.dbfprintf('\nDone finding events. Final total: %d\n', nQuakes);
- warning('on','MATLAB:table:ModifiedAndSavedVarnames');
- fetchOutputs(f);
- data = fileread(fn);
- delete(fn);
- if isvalid(myfig)
- close(myfig)
+
+ case "PayloadTooLarge"
+ disp(resp.Body.Data)
+
+ % as of 2018-12-14, INGV returns this result when limit is exceeded. these depend on the error message wording
+ nFound = double(extractBefore(extractAfter(resp.Body.Data,'the number of requested events is "'),'";'));
+ maxSearchLimit = double(extractBefore(extractAfter(resp.Body.Data,'your request must be less then "'),'" events.'));
+ if ~ismissing(maxSearchLimit)
+ warning("maximum number of events [%d] exceeded. atttempting to limit results", maxSearchLimit);
+ % try again, in chunks of maxSearchLimit
+ disp('* trying again while limiting results')
+ [resp, ok] = get_in_chunks(myuri, maxSearchLimit, nFound);
+ else
+ uOutput=[];
+ ok=false;
end
- else
- data = webread([baseurl 'query'], varargin{:},'format','text',options);
- end
- catch ME
- cancel(f)
- if isvalid(myfig)
- close(myfig)
- end
- switch ME.identifier
- %case 'MATLAB:webservices:CopyContentToDataStreamError'
- otherwise
- txt = 'An error occurred attempting to reach the FDSN web services';
- errordlg(sprintf('%s\n\n%s\n\nidentifier: ''%s''', txt, ME.message, ME.identifier),...
- 'Error retrieving data');
- end
- uOutput=[];
- ok=false;
+
+ otherwise
+ warning("there was some sort of problem. The response follows")
+ disp(resp)
+ warning(resp.Body.Data)
+ uOutput=[];
+ ok=false;
+ end
+ if ok
+ data = char(resp.Body.Data);
+ else
return
end
@@ -369,4 +302,59 @@
end
+function [resp,ok] = get_in_chunks(myuri, maxAllowed, expected)
+ % try again, in chunks of maxSearchLimit
+ starts=1:maxAllowed:expected;
+ uOutput='';
+ for n=starts
+ fprintf(' ** retrieving events %d to %d\n', n, min(n+maxAllowed-1, expected));
+ resp = get_low_level_fdsn_query(myuri + "&offset=" + n + "&limit="+ maxAllowed);
+ if resp.StatusCode == "OK"
+ if isempty(uOutput)
+ uOutput=string(resp.Body.Data);
+ else
+ uOutput=uOutput + newline + extractAfter(resp.Body.Data, newline);
+ end
+ ok=true;
+ else
+ ok=false;
+ break
+ end
+ end
+ if ok
+ resp.Body.Data=uOutput;
+ end
+end
+
+function resp = get_low_level_fdsn_query(uri)
+ U = matlab.net.URI(uri);
+ method = matlab.net.http.RequestMethod.GET;
+ type1 = matlab.net.http.MediaType('text/*');
+ acceptField = matlab.net.http.field.AcceptField(type1);
+ uafield=matlab.net.http.HeaderField('UserAgent','MATLAB 9.4.0.949201 (R2018a) Update 6 ZMAP/7.1');
+ contentTypeField = matlab.net.http.field.ContentTypeField('text/plain');
+ header = [acceptField contentTypeField uafield];
+ request = matlab.net.http.RequestMessage(method,header);
+
+ consumer=matlab.net.http.io.StringConsumer;
+ options=matlab.net.http.HTTPOptions(...'SavePayload',true ...
+ 'ProgressMonitorFcn',@MyProgressMonitorNew ...
+ ,'UseProgressMonitor',true ...
+ ,'ConnectTimeout',30 ... % in seconds
+ );
+ [resp,req,hist] = request.send(U,options,consumer);
+ %{
+ % if there is an error, it would be shown in hist.Response.Body.Data
+ ss=strsplit(string(resp.Body.Data),newline)';
+ numel(ss)
+
+ %%
+ f=fopen('junkk.dat','w');
+ fprintf(f,"%s",resp.Body.Data); %resp.Body.Payload
+ fclose(f);
+ %%
+ ZG.primeCatalog = import_fdsn_event(1,'junk.dat')
+ % ZmapMainWindow(ZG.primeCatalog)
+ %}
+end
diff --git a/sample_low_level_http_request.m b/sample_low_level_http_request.m
new file mode 100644
index 00000000..3e0a2bc8
--- /dev/null
+++ b/sample_low_level_http_request.m
@@ -0,0 +1,33 @@
+% dealing with fdsn data more-or-less directly
+% this is a temporary script, and should be rolled into the
+% FDSN import routines
+% with the anticipation that matlab.net.http.ProgressMonitor
+% will be used to provide feedback.
+
+% U = matlab.net.URI('http://service.iris.edu/fdsnws/event/1/query?starttime=2018-01-11T00:00:00&orderby=time&format=text&nodata=404');
+tic
+U = matlab.net.URI('http://service.iris.edu/fdsnws/event/1/query?starttime=2018-09-11T00:00:00&orderby=time&format=text&nodata=404');
+method = matlab.net.http.RequestMethod.GET;
+type1 = matlab.net.http.MediaType('text/*');
+acceptField = matlab.net.http.field.AcceptField([type1]);
+contentTypeField = matlab.net.http.field.ContentTypeField('text/plain');
+header = [acceptField contentTypeField];
+request = matlab.net.http.RequestMessage(method,header);
+
+consumer=matlab.net.http.io.StringConsumer;
+
+[resp,req,hist] = request.send(U,matlab.net.http.HTTPOptions('SavePayload',true,'ProgressMonitorFcn',@MyProgressMonitorNew,'UseProgressMonitor',true),consumer);
+% show(request)
+% resp.show
+
+% if there is an error, it would be shown in hist.Response.Body.Data
+ss=strsplit(string(char(resp.Body.Data')),newline)';
+numel(ss)
+%%
+f=fopen('junkk.dat','w');
+fprintf(f,"%s",resp.Body.Data); %resp.Body.Payload
+fclose(f);
+%%
+ZG.primeCatalog = import_fdsn_event(1,'junk.dat')
+% ZmapMainWindow(ZG.primeCatalog)
+toc
diff --git a/src/+XYfun/stressgrid.m b/src/+XYfun/stressgrid.m
index 37d3b254..b6e5f263 100644
--- a/src/+XYfun/stressgrid.m
+++ b/src/+XYfun/stressgrid.m
@@ -1,7 +1,7 @@
classdef stressgrid < ZmapHGridFunction
% STRESSGRID calculate stress grid for event that have Dip, DipDirection, and Rake
properties
-
+ calcmethod
end
properties(Constant)
diff --git a/src/@ZmapMainWindow/catalog_menu.m b/src/@ZmapMainWindow/catalog_menu.m
index ea3ce3be..296cf181 100644
--- a/src/@ZmapMainWindow/catalog_menu.m
+++ b/src/@ZmapMainWindow/catalog_menu.m
@@ -85,6 +85,7 @@ function cb_recall(~,~)
hh=msgbox_nobutton('The catalog has been recalled.','Recall Catalog');
hh.delay_for_close(1);
+ %obj.replot_all();
else
warndlg('No catalog is currently memorized','Recall Catalog');
end
diff --git a/src/@ZmapMainWindow/create_all_menus.m b/src/@ZmapMainWindow/create_all_menus.m
index 2eb7d112..0f5f10cc 100644
--- a/src/@ZmapMainWindow/create_all_menus.m
+++ b/src/@ZmapMainWindow/create_all_menus.m
@@ -455,7 +455,7 @@ function explore_catalog(~,~)
function create_decluster_menu(parent)
submenu = parent;% uimenu(parent,'Label','Decluster the catalog');
- uimenu(submenu,'Label','Decluster (Reasenberg)',MenuSelectedField(),@(~,~)ReasenbergDeclusterClass(obj.catalog),...
+ uimenu(submenu,'Label','Decluster (Reasenberg)',MenuSelectedField(),@cb_reasen,...
'Separator','on');
uimenu(submenu,'Label','Decluster (Gardner & Knopoff)',...
'Enable','off',...
@@ -465,4 +465,18 @@ function cb_declus_inp(~,~)
[out,nMethod]=declus_inp(obj.catalog)
error('declustered. now what to do with results?');
end
+
+ function cb_reasen(~,~)
+ rdc = ReasenbergDeclusterClass(obj.catalog, "CalcFinishedFcn",@update_window_with_declust);
+
+ function update_window_with_declust()
+ if ~isempty(rdc.declusteredCatalog)
+ msg.dbdisp('replacing the catalog')
+ obj.rawcatalog = rdc.declusteredCatalog;
+ obj.replot_all()
+ else
+ errordlg('Empty declustered catalog, Main window will not be updated')
+ end
+ end
+ end
end
\ No newline at end of file
diff --git a/src/@ZmapMainWindow/plot_base_events.m b/src/@ZmapMainWindow/plot_base_events.m
index 76d8d995..40d3c375 100644
--- a/src/@ZmapMainWindow/plot_base_events.m
+++ b/src/@ZmapMainWindow/plot_base_events.m
@@ -202,6 +202,12 @@ function cb_crop_to_selection(~, ~)
yl = [min(ol(:,2)) max(ol(:,2))];
obj.map_axes.XLim = xl;
obj.map_axes.YLim = yl;
+
+ % update the grid boundaries, too
+ obj.gridopt.AbsoluteGridLimits = [xl , yl];
+ % obj.Grid = obj.Grid.MaskWithShape(obj.shape);
+ obj.Grid = ZmapGrid(obj.Grid.Name, 'FromGridOptions',obj.gridopt,'Shape',obj.shape);
+
end
function cb_crop_to_axes(~, ~)
diff --git a/src/as_class/ZmapFunction.m b/src/as_class/ZmapFunction.m
index 67a076e1..31256d41 100644
--- a/src/as_class/ZmapFunction.m
+++ b/src/as_class/ZmapFunction.m
@@ -62,6 +62,7 @@
AutoShowPlots logical = true; % will plots be generated upon calculation?
InteractiveMode logical = true;
DelayProcessing logical = false;
+ CalcFinishedFcn function_handle = @do_nothing ; %call this once calculation ("do it") has finished.
end
properties(Constant,Abstract)
@@ -145,6 +146,7 @@
p.addParameter('AutoShowPlots', obj.AutoShowPlots)
p.addParameter('InteractiveMode', obj.InteractiveMode);
p.addParameter('DelayProcessing', obj.DelayProcessing);
+ p.addParameter('CalcFinishedFcn', obj.CalcFinishedFcn);
p.parse(varginstuff{:});
% assign values from parameter to this object
@@ -154,6 +156,7 @@
obj.AutoShowPlots = p.Results.AutoShowPlots;
obj.InteractiveMode = p.Results.InteractiveMode;
obj.DelayProcessing = p.Results.DelayProcessing;
+ obj.CalcFinishedFcn = p.Results.CalcFinishedFcn;
end
function StartProcess(obj)
diff --git a/src/cgr_utils/NDK.m b/src/cgr_utils/NDK.m
index 807ee6bb..446c2687 100644
--- a/src/cgr_utils/NDK.m
+++ b/src/cgr_utils/NDK.m
@@ -177,7 +177,7 @@
mt.Properties.VariableNames={'mrr', 'mtt', 'mff', 'mrt', 'mrf', 'mtf'};
c.MomentTensor= mt;
c.Dip=double(obj.allNDKs.Dip_NodalPlane1);
- c.DipDirection=double(obj.allNDKs.Strike_NodalPlane1)+ 90;
+ c.DipDirection = mod(double(obj.allNDKs.Strike_NodalPlane1)+ 90, 360 );
c.Rake=double(obj.allNDKs.Rake_NodalPlane1);
end
end
diff --git a/src/cgr_utils/gui/fdsn_chooser.mlapp b/src/cgr_utils/gui/fdsn_chooser.mlapp
index caf48e20..d7dd6796 100644
Binary files a/src/cgr_utils/gui/fdsn_chooser.mlapp and b/src/cgr_utils/gui/fdsn_chooser.mlapp differ
diff --git a/src/cgr_utils/place2coords.m b/src/cgr_utils/place2coords.m
new file mode 100644
index 00000000..3ea3a9df
--- /dev/null
+++ b/src/cgr_utils/place2coords.m
@@ -0,0 +1,36 @@
+function [res]=place2coords(city_and_country)
+ % ouput is either 'latlon' or 'box'
+ % box is returned as [Wlat, Elat, Slon, Nlon]
+ % uses openstreetmap to get values
+
+ queryparams = 'format=json';
+ % queryparams=[queryparams, '&polygon_geojson=1'];
+ if isempty(city_and_country)
+ error("need to provide city,country or city,state,country")
+ end
+ parts = strip(split(city_and_country,','));
+ queryparams = [queryparams, '&city=',parts{1}];
+
+ if numel(parts)==3
+ queryparams=[queryparams,'&state=',parts{2}];
+ parts{2}='';
+ end
+
+ if numel(parts{end})==2
+ queryparams=[queryparams,'&countrycode=',parts{end}];
+ else
+ queryparams=[queryparams,'&country=',parts{end}];
+ end
+
+ baseUrl = 'https://nominatim.openstreetmap.org/search/query?';
+
+ res = webread(['https://nominatim.openstreetmap.org/search/query?',queryparams]);
+ for j = 1 : numel(res)
+ bb = res(j).boundingbox;
+ res(j).boundingbox = struct('minLat',str2double(bb{1}),...
+ 'maxLat',str2double(bb{2}),...
+ 'minLon',str2double(bb{3}),...
+ 'maxLon',str2double(bb{4}));
+ end
+end
+
\ No newline at end of file
diff --git a/src/cgr_utils/plot_strike_dip.m b/src/cgr_utils/plot_strike_dip.m
new file mode 100644
index 00000000..552fa8df
--- /dev/null
+++ b/src/cgr_utils/plot_strike_dip.m
@@ -0,0 +1,69 @@
+function plot_strike_dip(catalog,ax)
+ % plot strike and dip for a catalog using DipDirection and Dip
+ % specify the axes
+
+ dipdir=catalog.DipDirection;
+ dip=catalog.Dip;
+ x0=catalog.Longitude;
+ y0=catalog.Latitude;
+
+ strikescale=0.5;
+ dipscale=0.3;
+ linecolor = [1 0 0];
+ limit_to_axes = true;
+
+ % figure out some scaling stuff
+ axun = get(ax,'Units');
+ set(ax,'Units','pixels')
+ p=fix(get(ax,'Position'));
+ set(ax,'Units',axun);
+
+ % TODO: do something with p and the xlim to come up with an appropriately sized symbol
+
+
+ delete(findobj(gcf,'Tag','diptest'))
+
+
+ strike = wrapTo360(dipdir-90);
+ scaling = ax.DataAspectRatio;
+ labels = cellstr( num2str(dip) );
+ labels=" "+labels;
+
+ is_horiz = dip==0;
+ is_vert = dip==90;
+
+ dx = cosd(strike) .* strikescale ./scaling(2);
+ dy = sind(strike) .* strikescale ./ scaling(1);
+ xx = ([x0,x0,x0] + [-dx, dx, nan(size(x0))])';
+ yy = ([y0,y0,y0] + [-dy, dy, nan(size(y0))])';
+ xx=xx(:);
+ yy=yy(:);
+
+ hold on;
+ plot(ax,xx,yy,'r','linewidth',1,'Tag','diptest','DisplayName','Strikes')
+
+ dipx=cosd(dipdir(:))* dipscale./scaling(2);
+ dipy=sind(dipdir(:))* dipscale./scaling(1);
+ xdx = ([x0,x0,x0] + [zeros(size(x0)) - is_vert .* dipx, dipx, nan(size(x0))])';
+ ydy = ([y0,y0,y0] + [zeros(size(y0)) - is_vert .* dipy, dipy, nan(size(y0))])';
+
+ xdx=xdx(:);
+ ydy=ydy(:);
+ plot(ax,xdx,ydy,'Color',linecolor,'linewidth',1.5,'Tag','diptest','DisplayName','Dips')
+ if any(is_horiz)
+ plot(ax,x0(is_horiz),y0(is_horiz),'o','Color',linecolor,'linewidth',1.5,'Tag','diptest','DisplayName','HorizDips')
+ end
+ if limit_to_axes
+ rangeidx=in_range(x0,xlim) & in_range(y0,ylim);
+ text(ax,x0(rangeidx), y0(rangeidx), labels(rangeidx), 'VerticalAlignment','top',...
+ 'HorizontalAlignment','left',...
+ 'Tag','diptest',...
+ 'Color',linecolor .* 0.66);
+
+ else
+ text(ax,x0, y0, labels, 'VerticalAlignment','top',...
+ 'HorizontalAlignment','left',...
+ 'Tag','diptest',...
+ 'Color',linecolor .* 0.66);
+ end
+end
\ No newline at end of file
diff --git a/src/pvals/MyPvalClass.m b/src/pvals/MyPvalClass.m
index 6b3a12d6..bacfb164 100644
--- a/src/pvals/MyPvalClass.m
+++ b/src/pvals/MyPvalClass.m
@@ -14,7 +14,7 @@
% % DOIT
% [P, P_stddev, C, C_stddev, K, K_stddev] = mypvc.mypval2m()
%
- %
+
properties
targetCerr = 0.001 % calculation stops once error in C is below this value
targetPerr = 0.001 % calculation stops once error in P is below this value
@@ -431,13 +431,11 @@ function pvalcat(catalog)
minDaysAfterMainshock = days(str2double(answer{2}));
valeg2 = str2double(answer{3});
- % cut catalog at mainshock time:
- l = catalog.Date > mainshockDate;
- catalog = catalog.subset(l); %keep events AFTER mainshock
+ % keep events AFTER mainshock
+ catalog = catalog.subset( catalog.Date > mainshockDate );
- % cat at selected magnitude threshold
- l = catalog.Magnitude >= obj.MinThreshMag;
- catalog = catalog.subset(l); %keep big-enough events
+ % keep big-enough events
+ catalog = catalog.subset( catalog.Magnitude >= obj.MinThreshMag );
ZG.hold_state2=true;
ctp=CumTimePlot(catalog);
@@ -465,7 +463,6 @@ function pvalcat(catalog)
tmin = min(timeSinceMainshock);
tmax = max(timeSinceMainshock);
- tint = [tmin tmax];
[pv, pstd, cv, cstd, kv, kstd, rja, rjb] = mypval2m(obj, eqDates, eqMags,'date',valeg2,CO);
@@ -524,6 +521,7 @@ function pvalcat(catalog)
ratac = numv ./ dursir;
+ tint = [tmin tmax];
frf = kv ./ ((tavg + cv).^pv);
frf2 = kv ./ ((days(tint) + cv).^pv);
diff --git a/src/thomas/decluster/ReasenbergDeclusterClass.m b/src/thomas/decluster/ReasenbergDeclusterClass.m
index 9c1501e4..e613244e 100644
--- a/src/thomas/decluster/ReasenbergDeclusterClass.m
+++ b/src/thomas/decluster/ReasenbergDeclusterClass.m
@@ -5,32 +5,46 @@
properties
taumin duration = days(1) % look ahead time for not clustered events
taumax duration = days(10) % maximum look ahead time for clustered events
- P = 0.95 % confidence level, that you are observing next event in sequence
+ P = 0.95 % confidence level that this is next event in sequence
xk = 0.5 % is the factor used in xmeff
- xmeff = 1.5 % effective" lower magnitude cutoff for catalog. it is raised by a factor xk*cmag1 during clusters
+
+ % effective lower magnitude cutoff for catalog.
+ % During clusteres, it is raised by a factor xk*cmag1
+ xmeff = 1.5
+
rfact = 10 % factor for interaction radius for dependent events
err = 1.5 % epicenter error
derr = 2 % depth error, km
%declustRoutine = "ReasenbergDeclus"
declusteredCatalog ZmapCatalog
replaceSequenceWithEquivMainshock logical = false
- clusterDetailsVariableName char = 'cluster_details' % if empty, clustering details will not be saved to workspace
- declusteredCatalogVariableName char = 'declustered_catalog' % if empty, catalog will not be saved to workspace
+
+ % if empty, clustering details will not be saved to workspace
+ clusterDetailsVariableName char = 'cluster_details'
+
+ % if empty, catalog will not be saved to workspace
+ declusteredCatalogVariableName char = 'declustered_catalog'
+ memorizeOriginalCatalog logical = true
end
properties(Constant)
PlotTag = "ReasenbergDecluster"
+
ParameterableProperties = ["taumin", "taumax", "P",...
"xk","xmeff","rfact","err","derr",..."declustRoutine",...
"replaceSequenceWithEquivMainshock",...
"clusterDetailsVariableName",...
- "declusteredCatalogVariableName"];
- References = 'Paul Reasenberg (1985) "Second -order Moment of Central California Seismicity", JGR, Vol 90, P. 5479-5495.';
+ "declusteredCatalogVariableName",...
+ "memorizeOriginalCatalog"];
+
+ References = ['Paul Reasenberg (1985) ',...
+ '"Second -order Moment of Central California Seismicity"',...
+ ', JGR, Vol 90, P. 5479-5495.'];
end
methods
- function obj=ReasenbergDeclusterClass(catalog, varargin)
+ function obj = ReasenbergDeclusterClass(catalog, varargin)
% BVALGRID
% obj = BVALGRID() takes catalog, grid, and eventselection from ZmapGlobal.Data
%
@@ -61,17 +75,31 @@ function InteractiveSetup(obj)
zdlg.AddEdit('rfact', 'Interation radius factor:', obj.rfact, 'RFACT>/b>factor for interaction radius for dependent events');
zdlg.AddEdit('err', 'Epicenter error', obj.err, 'Epicenter error');
zdlg.AddEdit('derr', 'Depth error', obj.derr, 'derrDepth error');
+ %zdlg.AddHeader('Output')
+ zdlg.AddCheckbox('replaceSequenceWithEquivMainshock','Replace clusters with equivalent event',...
+ obj.replaceSequenceWithEquivMainshock, {},...
+ 'Will replace each set of cluster earthquakes with a single event of equivalent Magnitude');
+ zdlg.AddEdit('clusterDetailsVariableName', 'Save Clusters to workspace as', ...
+ obj.clusterDetailsVariableName, 'if empty, then nothing will be separately saved');
+ zdlg.AddEdit('declusteredCatalogVariableName', 'Save Declustered catalog to workspace as', ...
+ obj.declusteredCatalogVariableName,'if empty, then nothing will be separately saved');
+ zdlg.AddCheckbox('memorizeOriginalCatalog', 'Memorize Original catalog after sucessful decluster:', ...
+ obj.memorizeOriginalCatalog, {}, 'Memorize original catalog prior to declustering');
+
+
- zdlg.Create('Name', 'Reasenberg Declustering','WriteToObj',obj,'OkFcn',@obj.declus);
+ zdlg.Create('Name', 'Reasenberg Declustering','WriteToObj',obj,'OkFcn',@obj.Calculate);
end
function Results = Calculate(obj)
calcFn = @obj.declus;
- [declustered_catalog, misc] = calcFn(obj);
+ [obj.declusteredCatalog, misc] = calcFn(obj);
if nargout == 1
- Results = declustered_catalog;
+ Results = obj.declusteredCatalog;
end
+
+ obj.CalcFinishedFcn();
end
@@ -121,14 +149,14 @@ function InteractiveSetup(obj)
%% calculate interaction_zone (1 value per event)
interactzone_main_km = 0.011*10.^(0.4* obj.RawCatalog.Magnitude); %interaction zone for mainshock
- interactzone_in_clust_km = obj.rfact * interactzone_main_km; %interaction zone if included in a cluster
+ interactzone_in_clust_km = obj.rfact * interactzone_main_km; %interaction zone if included in a cluster
tau_min = days(obj.taumin);
tau_max = days(obj.taumax);
%calculation of the eq-time relative to 1902
- eqtime=clustime(obj.RawCatalog);
+ eqtime = clustime(obj.RawCatalog);
%variable to store information whether earthquake is already clustered
clus = zeros(1,obj.RawCatalog.Count);
@@ -140,6 +168,7 @@ function InteractiveSetup(obj)
drawnow
declustering_start = tic;
%for every earthquake in catalog, main loop
+ confine_value = @(value, min_val, max_val) max( min( value, max_val), min_val);
for i = 1: (obj.RawCatalog.Count-1)
% "myXXX" refers to the XXX for this event
@@ -152,7 +181,7 @@ function InteractiveSetup(obj)
end
% variable needed for distance and timediff
- my_cluster=clus(i);
+ my_cluster = clus(i);
alreadyInCluster = my_cluster~=0;
% attach interaction time
@@ -160,37 +189,36 @@ function InteractiveSetup(obj)
if alreadyInCluster
if my_mag >= max_mag_in_cluster(my_cluster)
max_mag_in_cluster(my_cluster) = my_mag;
- idx_biggest_event_in_cluster(my_cluster)=i;
- look_ahead_days=tau_min;
+ idx_biggest_event_in_cluster(my_cluster) = i;
+ look_ahead_days = tau_min;
else
bgdiff = eqtime(i) - eqtime(idx_biggest_event_in_cluster(my_cluster));
look_ahead_days = clustLookAheadTime(obj.xk, max_mag_in_cluster(my_cluster), obj.xmeff, bgdiff, obj.P);
- look_ahead_days = min(look_ahead_days, tau_max);
- look_ahead_days = max(look_ahead_days, tau_min);
+ look_ahead_days = confine_value(look_ahead_days, tau_min, tau_max);
end
else
- look_ahead_days=tau_min;
+ look_ahead_days = tau_min;
end
%extract eqs that fit interation time window
- [~,ac]=timediff(i, look_ahead_days, clus, eqtime);
+ [~,ac] = timediff(i, look_ahead_days, clus, eqtime);
if ~isempty(ac) %if some eqs qualify for further examination
- rtest1=interactzone_in_clust_km(i);
- if look_ahead_days==obj.taumin
+ rtest1 = interactzone_in_clust_km(i);
+ if look_ahead_days == obj.taumin
rtest2 = 0;
else
- rtest2=interactzone_main_km(idx_biggest_event_in_cluster(my_cluster));
+ rtest2 = interactzone_main_km(idx_biggest_event_in_cluster(my_cluster));
end
- if alreadyInCluster % if i is already related with a cluster
- tm1 = clus(ac) ~= my_cluster; %eqs with a clustnumber different than i
+ if alreadyInCluster % if i is already related with a cluster
+ tm1 = clus(ac) ~= my_cluster; % eqs with a clustnumber different than i
if any(tm1)
- ac=ac(tm1);
+ ac = ac(tm1);
end
bg_ev_for_dist = idx_biggest_event_in_cluster(my_cluster);
else
@@ -203,22 +231,22 @@ function InteractiveSetup(obj)
%extract eqs that fit the spatial interaction time
sl0 = dist1<= rtest1 | dist2<= rtest2;
- if any(sl0) %if some eqs qualify for further examination
- ll=ac(sl0); %eqs that fit spatial and temporal criterion
- lla=ll(clus(ll)~=0); %eqs which are already related with a cluster
- llb=ll(clus(ll)==0); %eqs that are not already in a cluster
- if ~isempty(lla) %find smallest clustnumber in the case several
- sl1=min(clus(lla)); %numbers are possible
+ if any(sl0) %if some eqs qualify for further examination
+ ll = ac(sl0); %eqs that fit spatial and temporal criterion
+ lla = ll(clus(ll)~=0); %eqs which are already related with a cluster
+ llb = ll(clus(ll)==0); %eqs that are not already in a cluster
+ if ~isempty(lla) %find smallest clustnumber in the case several
+ sl1 = min(clus(lla)); %numbers are possible
if alreadyInCluster
- my_cluster= min([sl1,my_cluster]);
+ my_cluster = min([sl1, my_cluster]);
else
my_cluster = sl1;
end
if clus(i)==0
- clus(i)=my_cluster;
+ clus(i) = my_cluster;
end
- %merge all related clusters together in the cluster with the smallest number
- sl2 = lla(clus(lla)~=my_cluster);
+ % merge related clusters together into cluster with the smallest number
+ sl2 = lla(clus(lla) ~= my_cluster);
for j1 = [i,sl2]
if clus(j1) ~= my_cluster
clus(clus==clus(j1)) = my_cluster;
@@ -227,62 +255,63 @@ function InteractiveSetup(obj)
end
if my_cluster==0 %if there was neither an event in the interaction zone nor i, already related to cluster
- k=k+1; %
- my_cluster=k;
+ k = k+1; %
+ my_cluster = k;
clus(i) = my_cluster;
max_mag_in_cluster(my_cluster) = my_mag;
idx_biggest_event_in_cluster(my_cluster) = i;
end
- if size(llb)>0 %attach clustnumber to events not already related to a cluster
+ if size(llb)>0 % attach clustnumber to events yet unrelated to a cluster
clus(llb) = my_cluster * ones(1,length(llb)); %
end
- end %if ac
- end %if sl0
- end %for loop
+ end %if ac
+ end %if sl0
+ end
close(wai);
msg.dbfprintf('Declustering complete. It took %g seconds\n',toc(declustering_start));
%% this table contains all we need to know about the clusters. maybe.
details = table;
- details.Properties.UserData=struct;
- for j=1 : numel(obj.ParameterableProperties)
- details.Properties.UserData.(obj.ParameterableProperties(j))=obj.(obj.ParameterableProperties(j));
+ details.Properties.UserData = struct;
+ for j = 1 : numel(obj.ParameterableProperties)
+ details.Properties.UserData.(obj.ParameterableProperties(j)) = obj.(obj.ParameterableProperties(j));
end
- details.Properties.Description='Details for cluster as done by reasenberg declustering';
- details.eventNumber = (1:obj.RawCatalog.Count)';
- details.clusterNumber = clus(:);
- details.clusterNumber(details.clusterNumber==0)=missing;
- details.isBiggest=false(size(details.clusterNumber));
- details.isBiggest(idx_biggest_event_in_cluster)=true;
- details.Latitude = obj.RawCatalog.Latitude;
- details.Properties.VariableUnits(width(details))={'degrees'};
+ details.Properties.Description = 'Details for cluster, from reasenberg declustering';
+ details.eventNumber = (1:obj.RawCatalog.Count)';
+ details.clusterNumber = clus(:);
+ details.clusterNumber(details.clusterNumber==0) = missing;
+ details.isBiggest = false(size(details.clusterNumber));
+ details.isBiggest(idx_biggest_event_in_cluster) = true;
- details.Longitude = obj.RawCatalog.Longitude;
- details.Properties.VariableUnits(width(details))={'degrees'};
+ details.Latitude = obj.RawCatalog.Latitude;
+ details.Properties.VariableUnits(width(details)) = {'degrees'};
- details.Depth = obj.RawCatalog.Depth;
- details.Properties.VariableUnits(width(details))={'kilometers'};
+ details.Longitude = obj.RawCatalog.Longitude;
+ details.Properties.VariableUnits(width(details)) = {'degrees'};
- details.Magnitude = obj.RawCatalog.Magnitude;
+ details.Depth = obj.RawCatalog.Depth;
+ details.Properties.VariableUnits(width(details)) = {'kilometers'};
- details.MagnitudeType = obj.RawCatalog.MagnitudeType;
+ details.Magnitude = obj.RawCatalog.Magnitude;
- details.Date = obj.RawCatalog.Date;
+ details.MagnitudeType = obj.RawCatalog.MagnitudeType;
- details.InteractionZoneIfMain= interactzone_main_km;
- details.Properties.VariableUnits(width(details))={'kilometers'};
+ details.Date = obj.RawCatalog.Date;
+
+ details.InteractionZoneIfMain = interactzone_main_km;
+ details.Properties.VariableUnits(width(details)) = {'kilometers'};
details.InteractionZoneIfInClust = interactzone_in_clust_km;
- details.Properties.VariableUnits(width(details))={'kilometers'};
+ details.Properties.VariableUnits(width(details)) = {'kilometers'};
clusterFreeCatalog = obj.RawCatalog.subset(ismissing(details.clusterNumber));
%biggest_events_in_cluster = obj.RawCatalog.subset(details.isBiggest);
- outputcatalog=clusterFreeCatalog;
+ outputcatalog = clusterFreeCatalog;
if ~any(clus)
return
@@ -292,93 +321,66 @@ function InteractiveSetup(obj)
%build a matrix clust that stored clusters
[~, biggest_events_in_cluster, max_mag_in_cluster,~,~] = funBuildclu(obj.RawCatalog,idx_biggest_event_in_cluster,clus,max_mag_in_cluster);
- %% replace each cluster sequence with one event that summarizes it
+
+ % replace cluster sequences with summary events
if obj.replaceSequenceWithEquivMainshock
- ans_ = 'Replace';
+ equi = obj.equevent(details(~ismissing(details.clusterNumber),:)); % calculates equivalent events
+ if isempty(equi)
+ disp('No clusters in the catalog with this input parameters');
+ return;
+ end
+ tmpcat = cat(clusterFreeCatalog, equi); %new, unsorted catalog
+
else
- ans_ = 'No';
+ tmpcat = cat(clusterFreeCatalog, biggest_events_in_cluster); % builds catalog with biggest events instead
+ disp('Original mainshocks kept');
end
- if obj.InteractiveMode
- ans_ = questdlg('Replace mainshocks with equivalent events?',...
- 'Replace mainshocks with equivalent events?',...
- 'Replace','No',ans_ );
- end
-
- switch ans_
- case 'Replace'
- equi = obj.equevent(details(~ismissing(details.clusterNumber),:)); % calculates equivalent events
- if isempty(equi)
- disp('No clusters in the catalog with this input parameters');
- return;
- end
- tmpcat=cat(clusterFreeCatalog, equi); %new catalog, but not sorted
- case 'No'
- tmpcat=cat(clusterFreeCatalog, biggest_events_in_cluster); % builds catalog with biggest events instead
- disp('Original mainshocks kept');
-
- end
+ tmpcat.sort('Date');
+ tmpcat.Name = string(tmpcat.Name) + " (declust)";
- %% save the clustering details
- if obj.InteractiveMode
- answer=inputdlg('Save cluster information to workspace as:','Declustering details',...
- 1,{obj.clusterDetailsVariableName});
- else
- answer = {obj.clusterDetailsVariableName};
+ % save clustering details to workspace
+ if ~isempty(obj.clusterDetailsVariableName)
+ assignin('base',matlab.lang.makeValidName(obj.clusterDetailsVariableName),details);
end
- if ~isempty(answer)
- assignin('base',matlab.lang.makeValidName(answer{1}),details);
- end
-
- tmpcat.sort('Date')
-
- ZG = ZmapGlobal.Data;
- ZG.original=obj.RawCatalog; %save catalog in variable original
- ZG.newcat=ZG.primeCatalog;
- ZG.storedcat=ZG.original;
- ZG.cluscat=ZG.original.subset(clus(clus~=0));
-
- %% save the declustered catalog
- if obj.InteractiveMode
- answer=inputdlg('Save declustered catalog to workspace as:','Declustered catalog',...
- 1,{obj.declusteredCatalogVariableName});
- else
- answer = {obj.declusteredCatalogVariableName};
- end
- if ~isempty(answer)
- assignin('base', answer{1}, tmpcat);
+ if obj.memorizeOriginalCatalog
+ mm = MemorizedCatalogManager();
+ mm.memorize(obj.RawCatalog,'predeclust')
end
+ ZG = ZmapGlobal.Data;
+ ZG.original = obj.RawCatalog; %save catalog in variable original
+ %ZG.newcat = ZG.primeCatalog;
+ %ZG.storedcat = ZG.original;
+ ZG.cluscat = ZG.original.subset(clus(clus~=0));
- %{
- warning('should somehow zmap_update_displays()');
- plot_ax = findobj(gcf,'Tag','mainmap_ax');
- hold(plot_ax,'on');
- pl=scatter3(plot_ax,ZG.cluscat.Longitude, ZG.cluscat.Latitude,ZG.cluscat.Depth,[],'m', 'DisplayName','Clustered Events');
- pl.ZData=ZG.cluscat.Depth;
- %}
+ % save declustered catalog to workspace
+ if ~isempty(obj.declusteredCatalogVariableName)
+ assignin('base', obj.declusteredCatalogVariableName, tmpcat);
+ end
+
st1 = sprintf([' The declustering found %d clusters of earthquakes, a total of %d'...
- ' events (out of %d). The map window [would] now display the declustered catalog containing %d events.'...
- 'The individual clusters are displayed as magenta on the map.' ], ...
+ ' events (out of %d). The map window now displays the declustered catalog containing %d events.'], ...
biggest_events_in_cluster.Count, ZG.cluscat.Count, ZG.original.Count , ZG.primeCatalog.Count);
msgbox(st1,'Declustering Information')
watchoff
- outputcatalog=ZG.cluscat;
+ outputcatalog = tmpcat;
obj.Result(1).values.cluster_details = details;
+
end
function plot(obj, varargin)
- f=figure('Name','Reasenberg Deslustering Results');
- ax=subplot(2,2,1);
- ZG=ZmapGlobal.Data;
+ f = figure('Name','Reasenberg Deslustering Results');
+ ax = subplot(2,2,1);
+ ZG = ZmapGlobal.Data;
biggest = obj.Result.values.cluster_details(obj.Result.values.cluster_details.isBiggest,:);
- non_cluster= obj.Result.values.cluster_details(ismissing(obj.Result.values.cluster_details.clusterNumber),:);
- msf=str2func(ZG.MainEventOpts.MarkerSizeFcn);
+ non_cluster = obj.Result.values.cluster_details(ismissing(obj.Result.values.cluster_details.clusterNumber),:);
+ msf = str2func(ZG.MainEventOpts.MarkerSizeFcn);
scatter3(biggest.Longitude,biggest.Latitude,biggest.Depth,msf(biggest.Magnitude));
ax.ZDir='reverse';
title(ax,'Biggest events in each cluster');
@@ -392,7 +394,7 @@ function plot(obj, varargin)
ax = subplot(2,2,2);
- ax=subplot(2,1,2);
+ ax = subplot(2,1,2);
isInClust = ~ismissing(obj.Result.values.cluster_details.clusterNumber);
isNotBig = ~obj.Result.values.cluster_details.isBiggest;
@@ -401,7 +403,7 @@ function plot(obj, varargin)
ax.YDir='reverse';
hold on;
scatter(biggest.Date,biggest.Depth,msf(biggest.Magnitude),categorical(biggest.clusterNumber),'DisplayName','primary events in each cluster');
- cb=colorbar;
+ cb = colorbar;
cb.Label.String = 'Cluster #';
ax.YLabel.String = 'Depth [km]';
ax.XLabel.String = 'Date';
@@ -411,14 +413,14 @@ function plot(obj, varargin)
end
methods(Static)
- function h=AddMenuItem(parent,catalog)
+ function h = AddMenuItem(parent,catalog)
% create a menu item
label='Reasenberg Decluster';
- h=uimenu(parent,'Label',label,MenuSelectedField(), @(~,~)ReasenbergDeclusterClass(catalog));
+ h = uimenu(parent,'Label',label,MenuSelectedField(), @(~,~)ReasenbergDeclusterClass(catalog));
end
- function equi=equevent(tb)
+ function equi = equevent(tb)
% equevent calc equivalent event to cluster
% equi = equevent(catalog, cluster, bg)
% catalog : earthquake catalog
@@ -431,19 +433,19 @@ function plot(obj, varargin)
%
report_this_filefun();
- equi=ZmapCatalog;
+ equi = ZmapCatalog;
equi.Name='clusters';
if isempty(tb)
return
end
- j=0;
+ j = 0;
nClusts = max(tb.clusterNumber);
[elat, elon, edep, emag] = deal(nan(nClusts,1));
edate(nClusts,1) = datetime(missing);
- emagtype(nClusts,1)=categorical(missing);
+ emagtype(nClusts,1) = categorical(missing);
- for n=1 : max(tb.clusterNumber)
+ for n = 1 : max(tb.clusterNumber)
clust_events = tb(tb.clusterNumber==n,:);
if isempty(clust_events)
continue;
@@ -451,7 +453,7 @@ function plot(obj, varargin)
j = j + 1;
eqmoment = 10.^(clust_events.Magnitude .* 1.2);
- emoment=sum(eqmoment); %moment
+ emoment = sum(eqmoment); %moment
weights = eqmoment./emoment; %weightfactor
elat(j) = sum(clust_events.Latitude .* weights);
@@ -471,10 +473,11 @@ function plot(obj, varargin)
equi.Date = edate(1:j);
equi.Magnitude = emag(1:j);
equi.MagnitudeType = emagtype(1:j);
- equi.Depth=edep(1:j);
+ equi.Depth = edep(1:j);
[equi.Dip, equi.DipDirection, equi.Rake]=deal(nan(size(equi.Date)));
end
end
+
end