Skip to content

Commit

Permalink
Merge pull request #90 from kwschultz/master
Browse files Browse the repository at this point in the history 
Preparing for new release, doc changes
  • Loading branch information
@kwschultz
kwschultz committed Sep 14, 2015
2 parents a134f20 + b28740e commit 2b5aab4
Show file tree
Hide file tree
Showing 15 changed files with 821 additions and 222 deletions.
203 changes: 174 additions & 29 deletions PyVQ/pyvq/pyvq.py
View file

Large diffs are not rendered by default.

Binary file added doc/graphics/Event_0_plot.jpg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/graphics/Event_4406_crop_logo.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/graphics/VQ_Logo_fancy.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added doc/graphics/greens_fitting.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ...LCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ...LCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ..._ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ..._ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 1 addition & 1 deletion doc/vq.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ @ARTICLE{2005AGUFMNG21A..03Y
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

@article{Turcotte:2007up,
@article{Turcotte2007up,
author = {Turcotte, D and Holliday, J and Rundle, J},
title = {{BASS, an alternative to ETAS}},
journal = {Geophys. Res. Lett},
Expand Down
300 changes: 253 additions & 47 deletions doc/vq.tex
View file

Large diffs are not rendered by default.

35 changes: 18 additions & 17 deletions examples/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,8 @@ FOREACH(NPROC ${NUM_PROCS})

FOREACH(TAPER_IND RANGE ${NUM_TAPER})
LIST(GET TAPER_METHODS ${TAPER_IND} TAPER)
LIST(GET METHOD_MEANSLIP ${TAPER_IND} MEANSLIP)
LIST(GET METHOD_INTEREVENT ${TAPER_IND} INTEREVENT)
#LIST(GET METHOD_MEANSLIP ${TAPER_IND} MEANSLIP)
#LIST(GET METHOD_INTEREVENT ${TAPER_IND} INTEREVENT)
FOREACH(RES ${RESOLUTIONS})
SET(TEST_DIR ${CMAKE_CURRENT_BINARY_DIR}/PROCS${NPROC}/${TAPER}/)
FILE(MAKE_DIRECTORY ${TEST_DIR})
Expand Down Expand Up @@ -84,23 +84,24 @@ FOREACH(NPROC ${NUM_PROCS})
ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/")

# Check that the mean slip is near our expectations
IF (NOT ${MEANSLIP} EQUAL 0)
ADD_TEST(NAME test_slip_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR}
COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_slip ${MEANSLIP})
SET_TESTS_PROPERTIES (test_slip_P${NPROC}_${TAPER}_${RES} PROPERTIES
DEPENDS run_P${NPROC}_${TAPER}_${RES}
ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/")
ENDIF (NOT ${MEANSLIP} EQUAL 0)
# Schultz: For now removing until we come up with better tests
#IF (NOT ${MEANSLIP} EQUAL 0)
# ADD_TEST(NAME test_slip_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR}
# COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_slip ${MEANSLIP})
# SET_TESTS_PROPERTIES (test_slip_P${NPROC}_${TAPER}_${RES} PROPERTIES
# DEPENDS run_P${NPROC}_${TAPER}_${RES}
# ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/")
#ENDIF (NOT ${MEANSLIP} EQUAL 0)

# Check that the mean interevent time is near our expectations
IF (NOT ${INTEREVENT} EQUAL 0)

ADD_TEST(NAME test_interevent_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR}
COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_interevent ${INTEREVENT})
SET_TESTS_PROPERTIES (test_interevent_P${NPROC}_${TAPER}_${RES} PROPERTIES
DEPENDS run_P${NPROC}_${TAPER}_${RES}
ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/")
ENDIF(NOT ${INTEREVENT} EQUAL 0)
# Schultz: For now removing until we come up with better tests
#IF (NOT ${INTEREVENT} EQUAL 0)
# ADD_TEST(NAME test_interevent_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR}
# COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_interevent ${INTEREVENT})
# SET_TESTS_PROPERTIES (test_interevent_P${NPROC}_${TAPER}_${RES} PROPERTIES
# DEPENDS run_P${NPROC}_${TAPER}_${RES}
# ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/")
#ENDIF(NOT ${INTEREVENT} EQUAL 0)
ENDFOREACH(RES ${RESOLUTIONS})
ENDFOREACH(TAPER_IND RANGE ${NUM_TAPER})
ENDFOREACH(NPROC ${NUM_PROCS})
Expand Down
238 changes: 235 additions & 3 deletions quakelib/src/QuakeLibIO.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1901,6 +1901,21 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) {
out_file << "</Folder>\n";
out_file << "<Folder id=\"faults\">\n";


// Loop thru elements to compute min/max slip rates for color bounds
double max_rate = 0;
double min_rate = 0;
for (eit=_elements.begin(); eit!=_elements.end(); ++eit) {
max_rate = fmax( max_rate, eit->second.slip_rate());
min_rate = fmin( min_rate, eit->second.slip_rate());
}

/////// Bounds for color in blue to red range, white in the middle
double y_min = 0;
double y_max = 255;
double x_min = min_rate;
double x_max = max_rate;

// Go through the sections
for (fit=_sections.begin(); fit!=_sections.end(); ++fit) {
// And output the elements for each section
Expand All @@ -1927,16 +1942,41 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) {
lld[3] = lld[2];
lld[2] = c.convert2LatLon(xyz[2]+(xyz[1]-xyz[0]));
}

// Output the KML format polygon for this section

// Compute blue to red color (RGB)
// Keep red scale (min is white, max is red) so red=255
// Blue and green are equal and vary from (255 for min vals to 0 for max vals)
int blue, green;
int red = y_max;
if (eit->second.slip_rate() == 0) {
blue = 0;
green = 0;
red = 0;
} else {
int interp_color = (int) linear_interp(eit->second.slip_rate(), x_min, x_max, y_min, y_max);
blue = y_max - interp_color;
green = blue;
}

// Output the KML format polygon for this element
out_file << "\t\t<Placemark>\n";
out_file << "\t\t<description>\n";
out_file << "Fault name: " << fit->second.name() << "\n";
out_file << "Element #: " << eit->second.id() << "\n";
out_file << "Slip rate: " << c.m_per_sec2cm_per_yr(eit->second.slip_rate()) << " cm/year\n";
out_file << "Rake: " << c.rad2deg(eit->second.rake()) << " degrees\n";
out_file << "Aseismic: " << eit->second.aseismic() << "\n";
out_file << "\t\t</description>\n";
out_file << "\t\t\t<styleUrl>#baseStyle</styleUrl>\n";
out_file << "\t\t\t<Style>\n";
out_file << "\t\t\t\t<LineStyle>\n";
out_file << "\t\t\t\t\t<color>"<< rgb2hex(red, green, blue) <<"</color>\n";
out_file << "\t\t\t\t\t<width>1</width>\n";
out_file << "\t\t\t\t</LineStyle>\n";
out_file << "\t\t\t\t<PolyStyle>\n";
out_file << "\t\t\t\t\t<color>"<< rgb2hex(red, green, blue) <<"</color>\n";
out_file << "\t\t\t\t</PolyStyle>\n";
out_file << "\t\t\t</Style>\n";
//out_file << "\t\t\t<styleUrl>#baseStyle</styleUrl>\n";
out_file << "\t\t\t<Polygon>\n";
out_file << "\t\t\t\t<extrude>0</extrude>\n";
out_file << "\t\t\t\t<altitudeMode>relativeToGround</altitudeMode>\n";
Expand Down Expand Up @@ -1967,6 +2007,198 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) {
return 0;
}

// Schultz: Using these functions for color interpolation and RGB to HEX color conversions
double quakelib::ModelWorld::linear_interp(const double &x, const double &x_min, const double &x_max, const double &y_min, const double &y_max) const {
return ((y_max - y_min)/(x_max - x_min) * (x - x_min)) + y_min;
}

char* quakelib::ModelWorld::rgb2hex(const int r, const int g, const int b) const {
// Returning ABGR to work with KML format
char *buf;
size_t sz;
sz = snprintf(NULL, 0, "FF%02X%02X%02X", b,g,r);
buf = (char *)malloc(sz + 1);
snprintf(buf, sz+1, "FF%02X%02X%02X", b,g,r);
return buf;
}

int quakelib::ModelWorld::write_event_kml(const std::string &file_name, const quakelib::ModelEvent &event) {
std::ofstream out_file;
std::map<UIndex, ModelSection>::const_iterator fit;
std::map<UIndex, ModelElement>::const_iterator eit;
LatLonDepth min_bound, max_bound, center;
Vec<3> min_xyz, max_xyz;
double dx, dy, range, mag, mean_slip, ev_year;
unsigned int i, trigger, ev_num;
ElementIDSet involved_elements;
ElementIDSet::iterator it;

out_file.open(file_name.c_str());

get_bounds(min_bound, max_bound);
center = LatLonDepth(max_bound.lat() - (max_bound.lat()-min_bound.lat())/2,
max_bound.lon() - (max_bound.lon()-min_bound.lon())/2);
Conversion c(center);
min_xyz = c.convert2xyz(min_bound);
max_xyz = c.convert2xyz(max_bound);
dx = max_xyz[0]-min_xyz[0];
dy = max_xyz[1]-min_xyz[1];
range = fmax(dx, dy) * (1.0/tan(c.deg2rad(30)));

double max_depth = fabs(min_bound.altitude());

out_file << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
out_file << "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n";
out_file << "<Document>\n";
out_file << "<LookAt>\n";
out_file << "\t<latitude>" << center.lat() << "</latitude>\n";
out_file << "\t<longitude>" << center.lon() << "</longitude>\n";
out_file << "\t<altitude>0</altitude>\n";
out_file << "\t<range>" << range << "</range>\n";
out_file << "\t<tilt>0</tilt>\n";
out_file << "\t<heading>0</heading>\n";
out_file << "\t<altitudeMode>absolute</altitudeMode>\n";
out_file << "</LookAt>\n";
out_file << "<Style id=\"sectionLabel\">\n";
out_file << "\t<IconStyle>\n";
out_file << "\t\t<Icon>\n";
out_file << "\t\t\t<href>http://maps.google.com/mapfiles/kml/paddle/wht-blank.png</href>\n";
out_file << "\t\t</Icon>\n";
out_file << "\t</IconStyle>\n";
out_file << "</Style>\n";

//// Get relevant event info
involved_elements = event.getInvolvedElements();
trigger = event.getEventTrigger();
mean_slip = event.calcMeanSlip();
ev_year = event.getEventYear();
ev_num = event.getEventNumber();

// Write event info above triggering element
out_file << "\t<Placemark id=\"Event_" << ev_num << "_label\">\n";
out_file << "\t\t<name>" << "Event "<< ev_num << "</name>\n";
out_file << "\t\t<description>\n";
out_file << "Trigger: Element #"<< trigger << " \n";
out_file << "Mean slip [m] "<< mean_slip << " \n";
out_file << "Event year "<< ev_year << " \n";
out_file << "\t\t</description>\n";
out_file << "\t\t<styleUrl>#sectionLabel</styleUrl>\n";
out_file << "\t\t<Point>\n";
out_file << "\t\t\t<extrude>1</extrude>\n";
out_file << "\t\t\t<altitudeMode>relativeToGround</altitudeMode>\n";
// Find the deepest element for this section
UIndex best_vertex;
double min_altitude = DBL_MAX, cur_alt;
ModelVertex vert;

eit = _elements.find(trigger);
for (i=0; i<3; ++i) {
cur_alt = _vertices.find(eit->second.vertex(i))->second.lld().altitude();
if (cur_alt < min_altitude) {
min_altitude = cur_alt;
best_vertex = eit->second.vertex(i);
}
}

vert = _vertices.find(best_vertex)->second;
out_file << "\t\t\t<coordinates>" << vert.lld().lon() << "," << vert.lld().lat() << "," << max_depth + vert.lld().altitude() << "</coordinates>\n";
out_file << "\t\t</Point>\n";
out_file << "\t</Placemark>\n";


// Loop thru elements to compute min/max slip rates for color bounds
double max_slip = -DBL_MAX;
double min_slip = DBL_MAX;
for (it=involved_elements.begin(); it!=involved_elements.end(); ++it) {
max_slip = fmax( max_slip, event.getEventSlip(*it));
min_slip = fmin( min_slip, event.getEventSlip(*it));
}

/////// Bounds for color in blue to red range, white in the middle
double y_min = -255;
double y_max = 255;
double x_max = fmax(fabs(max_slip), fabs(min_slip));
double x_min = -x_max;

out_file << "<Folder id=\"elements\">\n";
// Go through the involved elements
for (it=involved_elements.begin(); it!=involved_elements.end(); ++it) {
LatLonDepth lld[4];
unsigned int i, npoints;

for (i=0; i<3; ++i) {
std::map<UIndex, ModelVertex>::const_iterator vit;
vit = _vertices.find( _elements.find(*it)->second.vertex(i) );
lld[i] = vit->second.lld();
}

// If this is a quad element, calculate the 4th implicit point
if (eit->second.is_quad()) {
Vec<3> xyz[3];

for (i=0; i<3; ++i) xyz[i] = c.convert2xyz(lld[i]);

lld[3] = lld[2];
lld[2] = c.convert2LatLon(xyz[2]+(xyz[1]-xyz[0]));
}

// Compute blue to red color (RGB)
int blue, green, red;
int interp_color = (int) linear_interp(event.getEventSlip(*it), x_min, x_max, y_min, y_max);
if (interp_color > 0) {
// More red
red = y_max;
blue = y_max - interp_color;
green = blue;
} else {
blue = y_max;
red = y_max - fabs(interp_color);
green = red;
}

// Output the KML format polygon for this element
out_file << "\t\t<Placemark>\n";
out_file << "\t\t<description>\n";
out_file << "Element #: " << *it << "\n";
out_file << "Event slip [m]: " << event.getEventSlip(*it) << "\n";
out_file << "\t\t</description>\n";
out_file << "\t\t\t<Style>\n";
out_file << "\t\t\t\t<LineStyle>\n";
out_file << "\t\t\t\t\t<color>"<< rgb2hex(red, green, blue) <<"</color>\n";
out_file << "\t\t\t\t\t<width>1</width>\n";
out_file << "\t\t\t\t</LineStyle>\n";
out_file << "\t\t\t\t<PolyStyle>\n";
out_file << "\t\t\t\t\t<color>"<< rgb2hex(red, green, blue) <<"</color>\n";
out_file << "\t\t\t\t</PolyStyle>\n";
out_file << "\t\t\t</Style>\n";
out_file << "\t\t\t<Polygon>\n";
out_file << "\t\t\t\t<extrude>0</extrude>\n";
out_file << "\t\t\t\t<altitudeMode>relativeToGround</altitudeMode>\n";
out_file << "\t\t\t\t<outerBoundaryIs>\n";
out_file << "\t\t\t\t\t<LinearRing>\n";
out_file << "\t\t\t\t\t\t<coordinates>\n";
npoints = (eit->second.is_quad() ? 4 : 3);

for (i=0; i<npoints+1; ++i) out_file << "\t\t\t\t\t\t\t" << lld[i%npoints].lon() << "," << lld[i%npoints].lat() << "," << max_depth + lld[i%npoints].altitude() << "\n";

out_file << "\t\t\t\t\t\t</coordinates>\n";
out_file << "\t\t\t\t\t</LinearRing>\n";
out_file << "\t\t\t\t</outerBoundaryIs>\n";
out_file << "\t\t\t</Polygon>\n";
out_file << "\t\t</Placemark>\n";

}

out_file << "</Folder>\n";
out_file << "</Document>\n";
out_file << "</kml>\n";

out_file.close();

return 0;
}


void quakelib::ModelSection::apply_remap(const ModelRemapping &remap) {
_data._id = remap.get_section_id(_data._id);
}
Expand Down
Loading

0 comments on commit 2b5aab4

Please sign in to comment.