Let @J in GMT/OGR files set -fg if implying degree coordinates#6057
Let @J in GMT/OGR files set -fg if implying degree coordinates#6057PaulWessel merged 6 commits intomasterfrom
Conversation
Basically, if a shapefile or GMT/OGR file has projection information that says it is a long,lat file then we should set -fg internally since we know it is in degrees.
src/gmt_io.c
Outdated
| GMT->current.io.ogr = GMT_OGR_FALSE; | ||
| return (false); | ||
| case 'e': k = 0; /* EPSG code */ | ||
| if (strstr (Jstring, "4326")) /* 4326 means we have lon/lat degree coordinates and thus should set -fig */ |
There was a problem hiding this comment.
Several other EPSGs are also lat/long
src/gmt_io.c
Outdated
| break; | ||
| case 'g': k = 1; break; /* GMT proj code. If given then data are projected */ | ||
| case 'p': k = 2; /* Proj.4 code */ | ||
| if (strstr (Jstring, "+proj=longlat")) /* +proj=longlat means we have lon/lat degree coordinates and thus should set -fig */ |
There was a problem hiding this comment.
"+pronj=latlong" is also valid
src/gmt_io.c
Outdated
| gmt_set_geographic (GMT, GMT_IN); | ||
| break; | ||
| case 'w': k = 3; /* OGR WKT representation */ | ||
| if (strstr (Jstring, "Degree")) /* UNIT as Degree means we have lon/lat degree coordinates and thus should set -fig */ |
There was a problem hiding this comment.
Nope. It would be degree but I think many many WKTs have UNIT["degree",0.0174532925199433.
AS I said, the only way I know out of this is to convert the WKT to PROJ4 and parse that.
|
Thanks, I will redo this a bit introducing some enums for the 4 cases and a separate function to make this decision. While at it, I see we initialize 3 of the 4 proj strings to degrees basically, but leave out the EPSG: Would it be OK to add (for
in that context? |
|
Dies this answer the question? |
|
The update places the decision after we have read all the 0-4 |
src/gmt_io.c
Outdated
| if (S->proj[GMTIO_GCS]) { /* See if we have degree units */ | ||
| if (strstr (S->proj[GMTIO_GCS], "UNIT[\"Degree\",0.0174532925199433]")) { | ||
| if (S->proj[GMTIO_WKT]) { /* See if we have degree units */ | ||
| if (strstr (S->proj[GMTIO_WKT], "UNIT[\"Degree\",0.0174532925199433]")) { |
There was a problem hiding this comment.
No. This is wrong. Se for example
julia> epsg2proj(2000)
"+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465 +units=m +no_defs"
julia> println(epsg2wkt(2000))
PROJCS["Anguilla 1957 / British West Indies Grid",
GEOGCS["Anguilla 1957",
DATUM["Anguilla_1957",
SPHEROID["Clarke 1880 (RGS)",6378249.145,293.465,
AUTHORITY["EPSG","7012"]],
AUTHORITY["EPSG","6600"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
There was a problem hiding this comment.
OK, I have removed that check. So currently no check on WKT but presumably any of the other 2-3 will be more helpful.
Basically, if a GMT/OGR file (or a shapefile converted to one under the hood) has projection information that says it is a lon,lat file then we can set -fg internally since we know it is in degrees.