Error while getting geographic coordinates of projected location

PyLith
21

The error indicates that the Proj library could not convert the coordinates of a point in the mesh coordinate system to geographic coordinates while initializing one of the materials. This suggests that there is an error in how you are specifying the coordinate system for your domain.

22

Does it mean that the projected x y z coordinate of my mesh is not right?

Alvina KK

23

The coordinates may be correct, but the specification of what that coordinate system is to PyLith/Proj may not be correct.

24

I converted my coordinate in QGIS using UTM47S and defined the coordinate system in all the cfg file with UTM47S.

Which part of the file should I check to fix this problem?

Alvina KK

25

How are you specifying the UTM47S projection in your .cfg files? Can you provide the coordinates of one of your points in longitude, latitude in WGS84 and in UTM47S so that we can verify the coordinate transformation is consistent with Proj?

26

Dear Brad,

This is how I assign the projection in my pylith.cfg :
[pylithapp.mesh_generator.reader]
coordsys = spatialdata.geocoords.CSGeoProj
coordsys.space_dim = 3
coordsys.datum_horiz = WGS84
coordsys.datum_vert = mean sea level
coordsys.projector.projection = utm
coordsys.projector.proj_options = +zone=47 +south

I also have attached my files here Error while getting geographic coordinates of projected location - #3 by alvinakkuncoro

This is one of my coordinate :
Lon,Lat = 92.05539,5.99347
x,y = 10247559.53417,663951.76772

Best Regards,
Alvina KK

27

The point 92.05539E, 5.99347N is in UTM zone 46N.

Proj.4 gives the following conversion if using UTM zone 47S:

echo "92.05539 5.99347" | proj +proj=utm +zone=47 +south
-270425.05	10667376.03

Inverting the projection gives back the origin point:

echo "-270425.05    10667376.03" | invproj +proj=utm +zone=47 +south -f "%12.4f"
     92.0554	      5.9935

Inverting the coordinates of your x and y yields:

echo " 10247559.53417 663951.76772" | invproj +proj=utm +zone=47 +south -f "%12.4f"
   -173.7178	    -24.1950

There appears to be something wrong with how you are projecting your coordinates.

28

Dear Brad,

I am sorry the previous coordinate I gave you was converted in mercator. So, it must be different.
My center coordinate is 99.5 , -3.851. It is in UTM zone 47S.
All my point is converted referred to this coordinate.
The converted coordinate in QGIS (utm47s) : 555513.381 , 9574298.167.
I have checked in proj (utm47s) : 555513.40 9574326.09

The previous coordinate (92.05539,5.99347),
Converted coordinate in QGIS (utm 47S) : -270424.6511 , 10667376
In proj (utm47s) : -270425.05 10667376.03

Best Regards,
Alvina KK

29

I downloaded your Exodus file. It has 4326 points where the coordinates are NaN. This appears to be causing the problem.

Note: It looks like your mesh is 1700 km x 1000 km x 200 km. If you use UTM coordinates, I would expect significant distortion because UTM was not designed for such large regions. I recommend using a more appropriate projection like Albers equal area.

30

Dear Brad,

Thank your for finding the problem.

May I know what can cause this NaN value after I created the mesh?

I have prepared the geometry in Mercator as well just to prepare this UTM problem can’t be solved.

Alvina K K

31

I have no idea where the NaNs in your mesh are coming from. I would try a much coarser mesh and see if the NaNs go away.

Here is the code you can use to check if there are NanS in the coordinates of your vertices:

#!/usr/bin/env nemesis

import numpy
import netCDF4

exodus = netCDF4.Dataset("mesh/mesh_tet_smooth.exo", 'r')
x = exodus.variables['coordx'][:]
y = exodus.variables['coordy'][:]
z = exodus.variables['coordz'][:]

xnan = numpy.isnan(x)
ynan = numpy.isnan(y)
znan = numpy.isnan(z)

print("Number of NaNs: x:%d, y:%d, z:%d" % (numpy.sum(xnan), numpy.sum(ynan), numpy.sum(znan)))