[cig-commits] [commit] devel, master: added smoothed etopo-4 topo/bathy model, and modified routines accordingly (88e40a4)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:00:26 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 88e40a4634a749a7420aa69db6d0224115df5a80
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Sat Jul 17 23:50:35 2004 +0000
added smoothed etopo-4 topo/bathy model, and modified routines accordingly
>---------------------------------------------------------------
88e40a4634a749a7420aa69db6d0224115df5a80
topo_bathy/check_topo_bathy_PPM_image.f90 | 82 ---------
...2ascii.f90 => convert_etopo2_bin2ascii_Sun.f90} | 6 +-
topo_bathy/filter_etopo_file.f90 | 86 ---------
topo_bathy/image_topo_bathy_no_smooth.jpg | Bin 0 -> 587385 bytes
topo_bathy/image_topo_bathy_smooth_10.jpg | Bin 0 -> 306143 bytes
topo_bathy/image_topo_bathy_smooth_7.jpg | Bin 0 -> 325214 bytes
topo_bathy/readme_etopo2.txt | 104 +++++++++++
topo_bathy/smooth_topo_bathy_PPM_image.f90 | 193 +++++++++++++++++++++
8 files changed, 300 insertions(+), 171 deletions(-)
diff --git a/topo_bathy/check_topo_bathy_PPM_image.f90 b/topo_bathy/check_topo_bathy_PPM_image.f90
deleted file mode 100644
index 0d6a5f8..0000000
--- a/topo_bathy/check_topo_bathy_PPM_image.f90
+++ /dev/null
@@ -1,82 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 3 . 2
-! --------------------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology July 2002
-!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
-!
-!=====================================================================
-
- program check_topo_bathy_PPM_image
-!
-!--- read topography and bathymetry ASCII file and create PPM image to check it
-!
-!--- compile and link with file topo_bathy.f90 from main SEM code
-!
- implicit none
-
- include "../../constants.h"
-
-! use integer array to store values
- integer ibathy_topo(NX_BATHY,NY_BATHY)
-
- integer grey(NX_BATHY,NY_BATHY)
- integer ivalue,icurrent_rec,ix,iz
-
-! read the topography file
- print *,'reading topo file'
- print *,'file used has a resolution of ',RESOLUTION_TOPO_FILE
- call read_topo_bathy_file(ibathy_topo)
-
- print *,'min and max of topography = ',minval(ibathy_topo),maxval(ibathy_topo)
-
-! creating image with grey levels
- grey = 255 * (ibathy_topo - minval(ibathy_topo)) / (maxval(ibathy_topo) - minval(ibathy_topo))
- where(grey < 1) grey = 1
- where(grey > 255) grey = 255
-
-! store image in PPM format with grey levels
-
-! creating the header
- open(unit=27,file='____tutu1',status='unknown')
- write(27,100)
- write(27,101)
- write(27,102) NX_BATHY,NY_BATHY
- write(27,103)
- close(27)
-
- 100 format('P5')
- 101 format('# creator DK')
- 102 format(i6,' ',i6)
- 103 format('255')
-
-! create the PPM image
- print *,'creating PPM image'
- open(unit=27,file='____tutu2',status='unknown',access='direct',recl=1)
-
- icurrent_rec = 1
-
- do iz = 1,NY_BATHY
- do ix = 1,NX_BATHY
-
-! write grey
- write(27,rec=icurrent_rec) achar(grey(ix,iz))
- icurrent_rec = icurrent_rec + 1
-
- enddo
- enddo
-
- close(27)
-
- call system('cat ____tutu1 ____tutu2 > image_topo_bathy.ppm ; rm -f ____tutu1 ____tutu2')
-
- end program check_topo_bathy_PPM_image
-
diff --git a/topo_bathy/convert_etopo2_bin2ascii.f90 b/topo_bathy/convert_etopo2_bin2ascii_Sun.f90
similarity index 94%
rename from topo_bathy/convert_etopo2_bin2ascii.f90
rename to topo_bathy/convert_etopo2_bin2ascii_Sun.f90
index a18dbb1..5971520 100644
--- a/topo_bathy/convert_etopo2_bin2ascii.f90
+++ b/topo_bathy/convert_etopo2_bin2ascii_Sun.f90
@@ -1,11 +1,11 @@
!=====================================================================
!
-! S p e c f e m 3 D G l o b e V e r s i o n 3 . 2
+! S p e c f e m 3 D G l o b e V e r s i o n 3 . 4
! --------------------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology July 2002
+! (c) California Institute of Technology August 2003
!
! A signed non-commercial agreement is required to use this program.
! Please check http://www.gps.caltech.edu/research/jtromp for details.
@@ -17,7 +17,7 @@
program convert_etopo2_bin2ascii
!
-!---- convert etopo-2 raw binary Sun file to ASCII - use on Sun GPS machines
+!---- convert etopo-2 raw binary Sun file to ASCII - use on Sun GPS machines at Caltech
!
implicit none
diff --git a/topo_bathy/filter_etopo_file.f90 b/topo_bathy/filter_etopo_file.f90
deleted file mode 100644
index f8176b2..0000000
--- a/topo_bathy/filter_etopo_file.f90
+++ /dev/null
@@ -1,86 +0,0 @@
-
- program filter_topo_file
-
-! smooth topo file using a window filter
-
- implicit none
-
- include "../../constants.h"
-
-! filter final surface using box filter
- integer, parameter :: SIZE_FILTER_ONE_SIDE = 2
-
-! impose min and max of topography and bathymetry
- integer, parameter :: MIN_TOPO = -5000
- integer, parameter :: MAX_TOPO = +5000
-
- integer ix,iy,ixconv,iyconv,i,ic,ixcount,iycount
- integer ixbox,iybox,ixval,iyval
-
- double precision rlon,rlat,rx,ry,a,b
- double precision sumval,value,dist,sigma
-
-! use integer array to store values
- integer itopo(NX_BATHY,NY_BATHY)
- integer itopo_filtered(NX_BATHY,NY_BATHY)
-
- integer itopo_x,itopo_y
-
- print *
- print *,'size of window filter is ',2*SIZE_FILTER_ONE_SIDE + 1
- print *
-
- print *,'imposing min and max of topography and bathymetry = ',MIN_TOPO,MAX_TOPO
- print *
-
- open(unit=13,file='topo_bathy_etopo5.dat',status='old')
- do itopo_y=1,NY_BATHY
- do itopo_x=1,NX_BATHY
- read(13,*) itopo(itopo_x,itopo_y)
- enddo
- enddo
- close(13)
-
- print *,'min and max of original topo file is ',minval(itopo),maxval(itopo)
-
-! impose min and max values
- do itopo_y=1,NY_BATHY
- do itopo_x=1,NX_BATHY
- if(itopo(itopo_x,itopo_y) < MIN_TOPO) itopo(itopo_x,itopo_y) = MIN_TOPO
- if(itopo(itopo_x,itopo_y) > MAX_TOPO) itopo(itopo_x,itopo_y) = MAX_TOPO
- enddo
- enddo
-
-! filter final surface using box filter
- print *,'filtering topography data file using box filter'
- do iy = 1,NY_BATHY
- print *,'doing iy = ',iy,' out of ',NY_BATHY
- do ix = 1,NX_BATHY
- sumval = 0.d0
- do iybox = iy-SIZE_FILTER_ONE_SIDE,iy+SIZE_FILTER_ONE_SIDE
- do ixbox = ix-SIZE_FILTER_ONE_SIDE,ix+SIZE_FILTER_ONE_SIDE
- ixval = ixbox
- iyval = iybox
- if(ixval < 1) ixval = NX_BATHY - abs(ixval)
- if(iyval < 1) iyval = NY_BATHY - abs(iyval)
- if(ixval > NX_BATHY) ixval = ixval - NX_BATHY
- if(iyval > NY_BATHY) iyval = iyval - NY_BATHY
- sumval = sumval + dble(itopo(ixval,iyval))
- enddo
- enddo
- itopo_filtered(ix,iy) = nint(sumval/dble((2*SIZE_FILTER_ONE_SIDE+1)**2))
- enddo
- enddo
-
- print *,'min and max of filtered topo file is ',minval(itopo_filtered),maxval(itopo_filtered)
-
- open(unit=13,file='topo_bathy_etopo5_filtered.dat',status='unknown')
- do itopo_y=1,NY_BATHY
- do itopo_x=1,NX_BATHY
- write(13,*) itopo_filtered(itopo_x,itopo_y)
- enddo
- enddo
- close(13)
-
- end program filter_topo_file
-
diff --git a/topo_bathy/image_topo_bathy_no_smooth.jpg b/topo_bathy/image_topo_bathy_no_smooth.jpg
new file mode 100644
index 0000000..e54c273
Binary files /dev/null and b/topo_bathy/image_topo_bathy_no_smooth.jpg differ
diff --git a/topo_bathy/image_topo_bathy_smooth_10.jpg b/topo_bathy/image_topo_bathy_smooth_10.jpg
new file mode 100644
index 0000000..042804b
Binary files /dev/null and b/topo_bathy/image_topo_bathy_smooth_10.jpg differ
diff --git a/topo_bathy/image_topo_bathy_smooth_7.jpg b/topo_bathy/image_topo_bathy_smooth_7.jpg
new file mode 100644
index 0000000..54e17b2
Binary files /dev/null and b/topo_bathy/image_topo_bathy_smooth_7.jpg differ
diff --git a/topo_bathy/readme_etopo2.txt b/topo_bathy/readme_etopo2.txt
new file mode 100644
index 0000000..3bc49f0
--- /dev/null
+++ b/topo_bathy/readme_etopo2.txt
@@ -0,0 +1,104 @@
+ETOPO2 Global 2-Minute Gridded Elevation Data Volume E1, 8/30/2001
+
+Description:
+ Topography for the entire world, including oceans, continents, and polar
+ regions at a 2 minute resolution (~4 km) in meters.
+
+GMT:
+ This data set is included in the grdraster dataset for ease of plotting
+
+File Format:
+ The raw data is in etopo2.raw
+ The file is a raw header-less binary stored as
+ 2-byte (16 bit) signed integers
+ 10800 columns and 5400 rows
+
+
+Raw to Gridded Data:
+ xyz2grd etopo2.raw -R-180/-90/-90/90 -I2m -Getopo2.grd -F -ZTLh -V
+ grd2xyz etopo2.grd -bo > etopo2.xyz
+
+
+Data Sources:
+ For information about where the data originated check out
+ et2src.htm and startet2.htm
+
+
+-----------------------------------------------------------------------
+
+
+This CD-ROM contains several elements: Software, Data, and Images:
+
+GEODAS (GEOphysical DAta System) is an interactive database management
+system developed by the National Geophysical Data Center (NGDC) for use
+in the assimilation, storage and retrieval of geophysical data. The
+GEODAS software is being used with several types of data including
+marine trackline geophysical data, hydrographic (bathymetric) survey
+data, aeromagnetic survey data, multibeam bathymetric data, and coastal
+relief model bathymetric gridded data.
+
+The GEODAS software which comes with this CD includes:
+ Browser access to grids -- open startet2.htm in your browser after
+ installing the GEODAS software from the "setup" directory
+ NGDC Grid Translator Program
+ Hydro-Plot Grid viewer Program
+
+Data:
+
+ Raw 2' grids in big-endian (ETOPO2.RAW) and little-endian (ETOPO2.dos)
+ formats
+ (16-bit signed integers, 10800 columns x 5400 rows)
+ See Major Sources of Data" on startet2.htm for details
+ For Geographic Information System (GIS) users, see the "ArcView Grids" link
+ in the GEODAS Help Index
+
+Images:
+
+ 4000x2000 color shaded relief JPEG image (grdet2.jpg)
+ and 45-degree shaded relief image tiles with 2' or 5' pixels (1350x1350 pixels
+ or 512x512 pixels), accessed through the STARTET2.HTM browser page (any browser
+ and computer type)
+
+
+********************
+*** INSTALLATION ***
+********************
+
+MS Windows 95/98/NT:
+To install see /setup/windows/readme.txt
+
+UNIX Xwindows:
+To install see /setup/xwindows/readme.txt
+
+Macintosh users must have a PC emulation program to use the GEODAS features;
+ the image display features of the STARTET2.HTM page work on any platform.
+
+****************************
+*** Technical Assistance ***
+****************************
+
+Dan Metzger
+Dan.R.Metzger at noaa.gov
+303-497-6542
+
+Peter Sloss
+Peter.W.Sloss at noaa.gov
+303-497-6119
+
+Dave Divins
+David.Divins at noaa.gov
+303-497-6505
+
+
+Get the latest version of the software at our web site:
+http://www.ngdc.noaa.gov/mgg/gdas/gx_announce.Html
+
+
+
+
+
+
+
+
+
+
diff --git a/topo_bathy/smooth_topo_bathy_PPM_image.f90 b/topo_bathy/smooth_topo_bathy_PPM_image.f90
new file mode 100644
index 0000000..38ea367
--- /dev/null
+++ b/topo_bathy/smooth_topo_bathy_PPM_image.f90
@@ -0,0 +1,193 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 3 . 4
+! --------------------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology August 2003
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+ program smooth_topo_bathy_PPM_image
+!
+!--- read topography and bathymetry ASCII file, smooth it and create PNM image to check it
+!
+
+!! DK DK this program uses a very simple window filter
+!! DK DK it works fine because the model has a very high resolution
+!! DK DK in principle it would be better to use a circular cap to take
+!! DK DK the notion of angular distance between points into account in the filter
+!! DK DK in practice though this simple filter works perfectly fine
+
+ implicit none
+
+ include "../../constants.h"
+
+! filter final surface using box filter
+ logical, parameter :: SMOOTH_THE_MODEL = .true.
+ integer, parameter :: SIZE_FILTER_ONE_SIDE = 7
+
+! use integer array to store values
+ integer ibathy_topo(NX_BATHY,NY_BATHY)
+ integer ibathy_topo_ori(NX_BATHY,NY_BATHY)
+
+ integer ivalue,icurrent_rec,ix,iy
+ integer minvalue,maxvalue
+ integer ix_current,iy_current,ix_min,ix_max,iy_min,iy_max,ix_value,iy_value
+
+ double precision value_sum,area_window
+
+!----
+
+! read the topography file
+ print *,'NX_BATHY,NY_BATHY = ',NX_BATHY,NY_BATHY
+ print *,'file used has a resolution of ',RESOLUTION_TOPO_FILE,' minutes'
+
+ print *
+ print *,'reading topo file'
+
+ open(unit=13,file='topo_bathy_etopo4_ori_from_etopo2_subsampled.dat',status='old')
+ do iy=1,NY_BATHY
+ do ix=1,NX_BATHY
+ read(13,*) ibathy_topo_ori(ix,iy)
+ enddo
+ enddo
+ close(13)
+
+!! DK DK compute min and max before smoothing
+ minvalue = minval(ibathy_topo_ori)
+ maxvalue = maxval(ibathy_topo_ori)
+ print *,'min and max of topography before smoothing = ',minvalue,maxvalue
+
+!----
+
+!! DK DK smooth topography/bathymetry model
+
+ if(SMOOTH_THE_MODEL) then
+
+ print *
+ print *,'smoothing topo file'
+ if(SIZE_FILTER_ONE_SIDE < 1) stop 'SIZE_FILTER_ONE_SIDE must be greater than 1 for filter'
+ print *,'size of window filter is ',2*SIZE_FILTER_ONE_SIDE+1,' x ',2*SIZE_FILTER_ONE_SIDE+1
+ area_window = dble((2*SIZE_FILTER_ONE_SIDE+1)**2)
+
+ do iy_current = 1,NY_BATHY
+
+ if(mod(iy_current,10) == 0) print *,'smoothing line ',iy_current,' out of ',NY_BATHY
+
+ do ix_current = 1,NX_BATHY
+
+ value_sum = 0.d0
+
+! compute min and max of window
+ ix_min = ix_current - SIZE_FILTER_ONE_SIDE
+ ix_max = ix_current + SIZE_FILTER_ONE_SIDE
+
+ iy_min = iy_current - SIZE_FILTER_ONE_SIDE
+ iy_max = iy_current + SIZE_FILTER_ONE_SIDE
+
+! loop on points in window to compute sum
+ do iy = iy_min,iy_max
+ do ix = ix_min,ix_max
+
+! copy current value
+ ix_value = ix
+ iy_value = iy
+
+! avoid edge effects, use periodic boundary
+ if(ix_value < 1) ix_value = ix_value + NX_BATHY
+ if(ix_value > NX_BATHY) ix_value = ix_value - NX_BATHY
+
+ if(iy_value < 1) iy_value = iy_value + NY_BATHY
+ if(iy_value > NY_BATHY) iy_value = iy_value - NY_BATHY
+
+! compute sum
+ value_sum = value_sum + dble(ibathy_topo_ori(ix_value,iy_value))
+
+ enddo
+ enddo
+
+! assign mean value to filtered point
+ ibathy_topo(ix_current,iy_current) = nint(value_sum / area_window)
+
+ enddo
+ enddo
+
+ else
+
+! no smoothing
+ ibathy_topo = ibathy_topo_ori
+
+ endif
+
+!----
+
+!! DK DK compute min and max after smoothing
+ minvalue = minval(ibathy_topo)
+ maxvalue = maxval(ibathy_topo)
+ print *,'min and max of topography after smoothing = ',minvalue,maxvalue
+
+!! DK DK save the smoothed model
+ if(SMOOTH_THE_MODEL) then
+ print *
+ print *,'saving the smoothed model'
+ open(unit=13,file='topo_bathy_etopo4_smoothed_window7.dat',status='unknown')
+ do iy=1,NY_BATHY
+ do ix=1,NX_BATHY
+ write(13,*) ibathy_topo(ix,iy)
+ enddo
+ enddo
+ close(13)
+ print *,'can suppress white spaces in filtered model to save space if needed'
+ endif
+
+! create image with grey levels
+ do iy = 1,NY_BATHY
+ do ix = 1,NX_BATHY
+ ibathy_topo(ix,iy) = 255 * (ibathy_topo(ix,iy) - minvalue) / (maxvalue - minvalue)
+ if(ibathy_topo(ix,iy) < 1) ibathy_topo(ix,iy) = 1
+ if(ibathy_topo(ix,iy) > 255) ibathy_topo(ix,iy) = 255
+ enddo
+ enddo
+
+! store image in PNM format with grey levels
+
+! create the PNM image
+ print *
+ print *,'creating PNM image'
+
+! creating the header
+ open(unit=27,file='image_topo_bathy.pnm',status='unknown')
+ write(27,100)
+ write(27,101)
+ write(27,102) NX_BATHY,NY_BATHY
+ write(27,103)
+
+ 100 format('P3')
+ 101 format('# creator DK')
+ 102 format(i6,' ',i6)
+ 103 format('255')
+ 104 format(i3)
+
+ do iy = 1,NY_BATHY
+ do ix = 1,NX_BATHY
+
+! write data value (red = green = blue to produce grey levels)
+ write(27,104) ibathy_topo(ix,iy)
+ write(27,104) ibathy_topo(ix,iy)
+ write(27,104) ibathy_topo(ix,iy)
+
+ enddo
+ enddo
+
+ close(27)
+
+ end program smooth_topo_bathy_PPM_image
+
More information about the CIG-COMMITS
mailing list