need more than UTM zone number (southern hemisphere)

I am trying a simulation for New Zealand -- my first attempt for a region in the southern hemisphere. Here is the results from output_solver.txt.

//////////////// original (requested) position of the source: latitude: -38.9200000000000 longitude: 178.400000000000 UTM x: 621370.527514857 UTM y: -4308615.74860358 depth: 30.4000000000000 km topo elevation: -2333.72827148438
position of the source that will be used: UTM x: 621370.527514857 UTM y: 5335486.00000000 depth: 40.8607170410156 km z: -43194.4453125000
error in location of the source: 9644107. m ////////////////

The UTM zone for this region is 60H. The Par_file only supports a number (60), and it appears to default to northern hemisphere. Here is a test using all the default ellipsoids in Matlab.

Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60H) --> 621349.69 5691557.05 bessel ellipsoid (utm zone 60H) --> 621352.86 5691600.63 airy ellipsoid (utm zone 60H) --> 621355.94 5691475.19 clarke66 ellipsoid (utm zone 60H) --> 621370.53 5691377.17 clarke80 ellipsoid (utm zone 60H) --> 621372.18 5691467.78 international ellipsoid (utm zone 60H) --> 621372.88 5691096.83 krasovsky ellipsoid (utm zone 60H) --> 621369.45 5691093.16 wgs60 ellipsoid (utm zone 60H) --> 621367.93 5691147.20 iau65 ellipsoid (utm zone 60H) --> 621367.86 5691154.42 wgs66 ellipsoid (utm zone 60H) --> 621367.58 5691164.55 iau68 ellipsoid (utm zone 60H) --> 621367.87 5691154.63 wgs72 ellipsoid (utm zone 60H) --> 621367.38 5691170.54 grs80 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 wgs84 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60S) --> 621349.69 -4308442.95 bessel ellipsoid (utm zone 60S) --> 621352.86 -4308399.37 airy ellipsoid (utm zone 60S) --> 621355.94 -4308524.81 clarke66 ellipsoid (utm zone 60S) --> 621370.53 -4308622.83 clarke80 ellipsoid (utm zone 60S) --> 621372.18 -4308532.22 international ellipsoid (utm zone 60S) --> 621372.88 -4308903.17 krasovsky ellipsoid (utm zone 60S) --> 621369.45 -4308906.84 wgs60 ellipsoid (utm zone 60S) --> 621367.93 -4308852.80 iau65 ellipsoid (utm zone 60S) --> 621367.86 -4308845.58 wgs66 ellipsoid (utm zone 60S) --> 621367.58 -4308835.45 iau68 ellipsoid (utm zone 60S) --> 621367.87 -4308845.37 wgs72 ellipsoid (utm zone 60S) --> 621367.38 -4308829.46 grs80 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60 wgs84 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60

You see that SPECFEM3D is returning values that are pretty close to the 60S values (S in the northern hemisphere), but way off from the correct 60H points (H in southern hemisphere).

BUG: UTM conversion will work for zone numbers that are in the northern hemisphere only. (Or maybe I missed some parameter in constants.h that allows the user to pick northern or southern hemisphere.) It might be safest to ask the user to list the utm zone in the Par_file as letter+number, like 11S or 60H.
QUESTION: What ellipsoid is being used? My numbers in the list above do not match the number from the SPECFEM3D output. The utmx matches that of clarke66, but the utmy does not match. It would be useful to check the UTM conversion with matlab or python "standard" formulas.

I have not opened the source code for this issue, but I'll look into more soon. At the moment, the simulation cannot run, since no receivers are (incorrectly) located within the region. Please chime in if you're familiar with the UTM conversion in the code. These details are not in the manual yet. (Thanks.)

Hi Carl, Hi all,

I took that routine (utm_geo.f90) from an open-source package called CAMx ( http://www.camx.com/ ) a long time ago (around 2001 / 2002). You are right, back then apparently there was a bug in their routine in the case of the Southern hemisphere because in their current version they have added this comment:

c Modifications: c 12/12/02 Added logic for southern hemisphere UTM zones c Use zones +1 to +60 for NH, -1 to -60 for SH c Equator is defined as 0 km North for NH, 10,000 km N for SH

thus let me replace the old one with the new one and we should be all set. I am unavailable today, I will do it tomorrow.

Thanks, Dimitri.

On 11/08/2014 00:17, Carl Tape wrote:

I am trying a simulation for New Zealand -- my first attempt for a region in the southern hemisphere. Here is the results from output_solver.txt.

original (requested) position of the source: latitude: -38.9200000000000

longitude: 178.400000000000

UTM x: 621370.527514857

UTM y: -4308615.74860358

depth: 30.4000000000000 km topo elevation: -2333.72827148438

position of the source that will be used: UTM x: 621370.527514857

UTM y: 5335486.00000000

depth: 40.8607170410156 km z: -43194.4453125000

error in location of the source: 9644107. m

The UTM zone for this region is 60H. The Par_file only supports a number (60), and it appears to default to northern hemisphere. Here is a test using all the default ellipsoids in Matlab.

Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60H) --> 621349.69 5691557.05 bessel ellipsoid (utm zone 60H) --> 621352.86 5691600.63 airy ellipsoid (utm zone 60H) --> 621355.94 5691475.19 clarke66 ellipsoid (utm zone 60H) --> 621370.53 5691377.17 clarke80 ellipsoid (utm zone 60H) --> 621372.18 5691467.78 international ellipsoid (utm zone 60H) --> 621372.88 5691096.83 krasovsky ellipsoid (utm zone 60H) --> 621369.45 5691093.16 wgs60 ellipsoid (utm zone 60H) --> 621367.93 5691147.20 iau65 ellipsoid (utm zone 60H) --> 621367.86 5691154.42 wgs66 ellipsoid (utm zone 60H) --> 621367.58 5691164.55 iau68 ellipsoid (utm zone 60H) --> 621367.87 5691154.63 wgs72 ellipsoid (utm zone 60H) --> 621367.38 5691170.54 grs80 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 wgs84 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60S) --> 621349.69 -4308442.95 bessel ellipsoid (utm zone 60S) --> 621352.86 -4308399.37 airy ellipsoid (utm zone 60S) --> 621355.94 -4308524.81 clarke66 ellipsoid (utm zone 60S) --> 621370.53 -4308622.83 clarke80 ellipsoid (utm zone 60S) --> 621372.18 -4308532.22 international ellipsoid (utm zone 60S) --> 621372.88 -4308903.17 krasovsky ellipsoid (utm zone 60S) --> 621369.45 -4308906.84 wgs60 ellipsoid (utm zone 60S) --> 621367.93 -4308852.80 iau65 ellipsoid (utm zone 60S) --> 621367.86 -4308845.58 wgs66 ellipsoid (utm zone 60S) --> 621367.58 -4308835.45 iau68 ellipsoid (utm zone 60S) --> 621367.87 -4308845.37 wgs72 ellipsoid (utm zone 60S) --> 621367.38 -4308829.46 grs80 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60 wgs84 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60

You see that SPECFEM3D is returning values that are pretty close to the 60S values (S in the northern hemisphere), but way off from the correct 60H points (H in southern hemisphere).

1.
BUG: UTM conversion will work for zone numbers that are in the
northern hemisphere only. (Or maybe I missed some parameter in
constants.h that allows the user to pick northern or southern
hemisphere.) It might be safest to ask the user to list the utm zone
in the Par_file as letter+number, like 11S or 60H.
2.
QUESTION: What ellipsoid is being used? My numbers in the list above
do not match the number from the SPECFEM3D output. The utmx matches
that of clarke66, but the utmy does not match. It would be useful to
check the UTM conversion with matlab or python "standard" formulas.
I have not opened the source code for this issue, but I'll look into more soon. At the moment https://github.com/geodynamics/specfem3d/issues/new# the simulation cannot run, since no receivers are (incorrectly) located within the region. Please chime in if you're familiar with the UTM conversion in the code. These details are not in the manual yet. (Thanks.)

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/212.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

Hi Carl,

If you want more details about how they implemented it, here is what they say:

! UTM_GEO performs UTM to geodetic (long/lat) translation, and back. ! ! This is a Fortran version of the BASIC program "Transverse Mercator ! Conversion", Copyright 1986, Norman J. Berls (Stefan Musarra, 2/94) ! Based on algorithm taken from "Map Projections Used by the USGS" ! by John P. Snyder, Geological Survey Bulletin 1532, USDI.

Best wishes, Dimitri.

On 11/08/2014 12:49, Dimitri Komatitsch wrote:

Hi Carl, Hi all,

I took that routine (utm_geo.f90) from an open-source package called CAMx ( http://www.camx.com/ ) a long time ago (around 2001 / 2002). You are right, back then apparently there was a bug in their routine in the case of the Southern hemisphere because in their current version they have added this comment:

c Modifications: c 12/12/02 Added logic for southern hemisphere UTM zones c Use zones +1 to +60 for NH, -1 to -60 for SH c Equator is defined as 0 km North for NH, 10,000 km N for SH

thus let me replace the old one with the new one and we should be all set. I am unavailable today, I will do it tomorrow.

Thanks, Dimitri.

On 11/08/2014 00:17, Carl Tape wrote:
I am trying a simulation for New Zealand -- my first attempt for a region in the southern hemisphere. Here is the results from output_solver.txt.

original (requested) position of the source: latitude: -38.9200000000000

longitude: 178.400000000000

UTM x: 621370.527514857

UTM y: -4308615.74860358

depth: 30.4000000000000 km topo elevation: -2333.72827148438

position of the source that will be used: UTM x: 621370.527514857

UTM y: 5335486.00000000

depth: 40.8607170410156 km z: -43194.4453125000

error in location of the source: 9644107. m

The UTM zone for this region is 60H. The Par_file only supports a number (60), and it appears to default to northern hemisphere. Here is a test using all the default ellipsoids in Matlab.

Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60H) --> 621349.69 5691557.05 bessel ellipsoid (utm zone 60H) --> 621352.86 5691600.63 airy ellipsoid (utm zone 60H) --> 621355.94 5691475.19 clarke66 ellipsoid (utm zone 60H) --> 621370.53 5691377.17 clarke80 ellipsoid (utm zone 60H) --> 621372.18 5691467.78 international ellipsoid (utm zone 60H) --> 621372.88 5691096.83 krasovsky ellipsoid (utm zone 60H) --> 621369.45 5691093.16 wgs60 ellipsoid (utm zone 60H) --> 621367.93 5691147.20 iau65 ellipsoid (utm zone 60H) --> 621367.86 5691154.42 wgs66 ellipsoid (utm zone 60H) --> 621367.58 5691164.55 iau68 ellipsoid (utm zone 60H) --> 621367.87 5691154.63 wgs72 ellipsoid (utm zone 60H) --> 621367.38 5691170.54 grs80 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 wgs84 ellipsoid (utm zone 60H) --> 621367.42 5691169.40 Point 1/1 at (178.400000, -38.920000) everest ellipsoid (utm zone 60S) --> 621349.69 -4308442.95 bessel ellipsoid (utm zone 60S) --> 621352.86 -4308399.37 airy ellipsoid (utm zone 60S) --> 621355.94 -4308524.81 clarke66 ellipsoid (utm zone 60S) --> 621370.53 -4308622.83 clarke80 ellipsoid (utm zone 60S) --> 621372.18 -4308532.22 international ellipsoid (utm zone 60S) --> 621372.88 -4308903.17 krasovsky ellipsoid (utm zone 60S) --> 621369.45 -4308906.84 wgs60 ellipsoid (utm zone 60S) --> 621367.93 -4308852.80 iau65 ellipsoid (utm zone 60S) --> 621367.86 -4308845.58 wgs66 ellipsoid (utm zone 60S) --> 621367.58 -4308835.45 iau68 ellipsoid (utm zone 60S) --> 621367.87 -4308845.37 wgs72 ellipsoid (utm zone 60S) --> 621367.38 -4308829.46 grs80 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60 wgs84 ellipsoid (utm zone 60S) --> 621367.42 -4308830.60

You see that SPECFEM3D is returning values that are pretty close to the 60S values (S in the northern hemisphere), but way off from the correct 60H points (H in southern hemisphere).

1.
BUG: UTM conversion will work for zone numbers that are in the
northern hemisphere only. (Or maybe I missed some parameter in
constants.h that allows the user to pick northern or southern
hemisphere.) It might be safest to ask the user to list the utm zone
in the Par_file as letter+number, like 11S or 60H.
2.
QUESTION: What ellipsoid is being used? My numbers in the list above
do not match the number from the SPECFEM3D output. The utmx matches
that of clarke66, but the utmy does not match. It would be useful to
check the UTM conversion with matlab or python "standard" formulas.
I have not opened the source code for this issue, but I'll look into more soon. At the moment https://github.com/geodynamics/specfem3d/issues/new# the simulation cannot run, since no receivers are (incorrectly) located within the region. Please chime in if you're familiar with the UTM conversion in the code. These details are not in the manual yet. (Thanks.)

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/212.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

Hi Dimitri,

I checked a couple things. utm_geo.f90 is indeed using Clarke 1866, as suggested from my test before. A check from Matlab below gives the same semimajor and semiminor axes listed in utm_geo.f90 (L71): double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0

referenceEllipsoid('clarke66','km') Properties: Code: 7008 Name: 'Clarke 1866' LengthUnit: 'kilometer' SemimajorAxis: 6378.2064 SemiminorAxis: 6356.5838 InverseFlattening: 294.978698213898 Eccentricity: 0.0822718542230038

As I mentioned before, though, for my test point the utmx value is exactly the same (matlab vs specfem), but the utmy value differs very slightly.

utmx_spec = 621370.527514857; utmy_spec = -4308615.74860358; utmx_matlab = 621370.527515; utmy_matlab = -4308622.832969; utmx_spec / utmx_matlab = 1.000000000000 utmy_spec / utmy_matlab = 0.999998355770

This is probably small enough to not matter for what we are doing, but it's worth listing. (The difference is 7 meters for the y coordinate.) I'm just checking the matlab mapping toolbox vs what we are using, so there's no guarantee that one or the other is correct.

The main message is that if users are making their own STATIONS files in utm and then converting to lon-lat, it is best to use the clarke66 ellipsoid to be consistent with SPECFEM3D. To change the ellipsoid used in the utm conversion in SPECFEM3D, one would need to change semimaj and semimin in utm_geo.f90.

Thanks, Carl

Switched to the new version of src/shared/utm_geo.f90 in order to have support for the Southern hemisphere. Also switched to the WGS84 ellipsoid instead of Clarke 1866. Updated the users manual accordingly.

60H

I only just read up on this today, but technically this is not UTM; it's MGRS.

The main fix to add the option of a negative sign in front of the UTM zone works. Here is a check for the updated code, which uses WGS84. This is for the input lon 178.40 lat -38.92, as before (utm zone -60).

output from specfem3d--> 621367.422820 5691176.332066 output from matlab [wgs84 ellipsoid, utm zone 60H] --> 621367.422820 5691169.402792 utmx_matlab / utmx_spec = 1.000000000000 utmy_matlab / utmy_spec = 0.999998782453

There are still some very small differences for the y-coordinate (about 7 meters).

SPECFEM / specfem3d

need more than UTM zone number (southern hemisphere) #212