UTM over different time zones and shape files
Closed this issue · 8 comments
Also known as the Ol Doinyo conundrum (the volcano sits between 36 and 37 South): UTM data cannot tackle a change in time zones - for example passing from 32 to 33 in the Alps. Without this, the only way is to go back to Cartesian coordinates. This sums up the issue with shapefiles. While Shaperead.jl can read them but there is no easy way to write them down. To talk with geomorphologists, would be a great upgrade to link to GDal.jl, which (to me) results quite complicated.
For now, my advice is to stick to cartesian if you want to include shapefiles or pass to lat-long if you go over different time zones. I can work on this issue but I do not see any easy fix.
If I am wrong, it would be good to add a tutorial on how to use UTM with different time zones.
I also had a look at ArchGDAL.jl, but adding it already produced a load of warnings, so I am not sure how well this works. Maybe we should resort for now to being able to read shape files, but not to write them, as you suggested.
Well the Oldoinyo Volcano can now be downloaded by specifying Lat/Lon coordinates (note that you have to use negative numbers for the southern hemisphere).
julia> Topo= ImportTopo([35.8, 36, -2.9, -2.7],file="@earth_relief_01s")
GeoData
size : (721, 721, 1)
lon ϵ [ 35.8 : 36.0]
lat ϵ [ -2.9 : -2.7]
depth ϵ [ 0.733 km : 3.238 km]
fields: (:Topography,)
julia> Write_Paraview(Topo,"Topo")
1-element Vector{String}:
"Topo.vts"which gives the following picture:
Next, you can convert the data to UTM using:
julia> Topo_UTM = convert(UTMData,Topo)
UTMData
UTM zone : 36-37 South
size : (721, 721, 1)
EW ϵ [ 166390.40523059748 : 833578.6766772219]
NS ϵ [ 9.679018014185755e6 : 9.701208067904716e6]
depth ϵ [ 733.0 m : 3238.0 m]
fields : (:Topography,)
which indeed shows that is in both UTM zones 35 and 36.
When it comes to shapefiles, I looked a bit at the ones you use in your volcano tutorial. They do have the UTM data of the coordinates, and seem to consist of a lot of polygons if loaded with Shapefile.jl. In principle, we should be able to convert that to Lon/Lat, but it seems to consist of deep layers of vectors. I have the impression that these are individual polygons, so once we extract that, it should be possible to save this as polydata using WriteVTK.jl.
I am still unable to install GMT.jl on my mac, despite having followed the installation procedures. Could you send me the topography between [35.75 36.3] and [-2.96 -2.44]? It is for Miriam's paper.
Regarding shape file, yes, I have also loaded them in Matlab with "shaperead("Morfo_Stru_CF_2016.shp"); - I attach a screenshot. They are poligons and we can save it with a different format for now.
Adding to this, the procedure is not straightforward for points, like earthquakes, see this:
julia> data = readdlm("locations.txt", ',', Float64, skipstart=0, header=false);
julia> lat = data[:,1];
julia> lon = data[:,2];
julia> depth = -data[:,3]/1000;
julia> EQ_Data_Degrees = GeophysicalModelGenerator.GeoData(lon, lat, depth, (Depth = depth * km,));
julia> Write_Paraview(EQ_Data_Degrees, "Ol_Earthquakes_Degrees", PointsData=true)
julia> EQ_Data_UTM = convert(UTMData,EQ_Data_Degrees)
julia> Write_Paraview(EQ_Data_UTM, "Ol_Earthquakes_UTM", PointsData=true)The result is this:
The earthquakes are divided between zones 36 (coordinates in the 80000s) and 37 (16000s) - any easy fix? I can obviously rescale manually, just wanted to know if there is a way to get around it.
No problem, Here the MWE:
julia> data = [-2.65133 36.148 -7.65649; -2.74147 35.9833 -10.5349];
julia> lon = data[:,1];
julia> lat = data[:,2];
julia> depth = data[:,3];
julia> EQ_Data_Degrees = GeophysicalModelGenerator.GeoData(lon, lat, depth, (Depth = depth * km,));
julia> Write_Paraview(EQ_Data_Degrees, "Ol_Earthquakes_Degrees", PointsData=true)
julia> EQ_Data_UTM = convert(UTMData,EQ_Data_Degrees);
julia> Write_Paraview(EQ_Data_UTM, "Ol_Earthquakes_UTM", PointsData=true)which looks like this.
h
Actually, I now realize that it is not a bug but in fact expected behaviour. If you translate a GeoData set to UTMData it will use the local UTM zone of every data point. In this case, they are in different zones, so the values look weird.
To solve this, you have two options, both of which require you to specify a ProjectionPoint, which is the point that you fix for your projection (from one map to the other).
julia> proj = ProjectionPoint(Lat=-2.7, Lon=36)Note: in the example above, I think you mixed up lat and lon, which confused me initially.
Should be
julia> data = [-2.65133 36.148 -7.65649; -2.74147 35.9833 -10.5349];
julia> lat = data[:,1];
julia> lon = data[:,2];
julia> depth = data[:,3];
julia> EQ_Data_Degrees = GeophysicalModelGenerator.GeoData(lon, lat, depth, (Depth = depth * km,));- Use only a single UTM zone for the projection, which can be done with
julia> EQ_Data_UTM = Convert2UTMzone(EQ_Data_Degrees, proj)
UTMData
UTM zone : 37-37 South
size : (2,)
EW ϵ [ 164543.07782450126 : 182849.94109182025]
NS ϵ [ 9.696560188501557e6 : 9.706580913421948e6]
depth ϵ [ -10534.9 m : -7656.49 m]
fields : (:Depth,)- Convert it to a
CartDatastructure that has km's in the 3 directions, using
julia> EQ_Data_CartData = Convert2CartData(EQ_Data_Degrees, proj)
CartData
size : (2,)
x ϵ [ -1.847327406096214 km : 16.459535861222772 km]
y ϵ [ -4.594761423205957 km : 5.425963497184217 km]
z ϵ [ -10.5349 km : -7.65649 km]
fields : (:Depth,)The ProjectionPoint will be (0,0) in the CartData data set.
I'm closing this as it is working as expected