ArctosDB / arctos

Arctos is a museum collections management system
https://arctos.database.museum
Apache License 2.0
59 stars 13 forks source link

Arctos Coordinate Converter #5910

Closed dustymc closed 1 year ago

dustymc commented 1 year ago

Is your feature request related to a problem? Please describe.

There are some scattered issues involving the coordinate converter. We can't always find or properly understand them, so consolidating here. This issue should relate only to the conversion process; we can sort out the rest once that's solidified.

Describe what you're trying to accomplish

Debug the Arctos Coordinate Converter

Describe the solution you'd like

Happy spatial data

Describe alternatives you've considered

HOW is this not a common thing?? Why am I writing this code and not just reusing something that's been in production for decades? Nobody seems to know, but if someone has code that'll take various formats (including UTM) and datums (horizontal only - easy peasy!) and magic it into standard coordinates then we can just use that.

Grittybits

The converter is https://github.com/ArctosDB/PG_DDL/blob/master/function/convertRawCoords.sql

It is essentially two operations (now documented with giant unmissable comments)

I've added a bunch of debug code, which hopefully will make it possible for someone to spot what's happening.

I've got this set up so that I can easily add more debug code (and hopefully use the function without changing that, but this is at test at the moment), so please let me know if some mystery remains and I'll try to fill in the gaps.

I have not figured out how to make this visible, except from a terminal, so for now I'll just run whatever's requested. Up for better ideas, and of course I'd welcome more users with sql access.

Please let me know if you want me to run some specific data, have a specific example to explore, whatever.

PLEASE let me know if you see anything in the examples and tests that doesn't make sense.

I'll use another comment for the first example, hopefully that'll let me keep the overview from being overwhelmed by the details.

@mkoo

dustymc commented 1 year ago

First test: UTM

Call


select convertRawCoords('{
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"13S",
    "utm_ns":"3885437",
    "utm_ew":"636394",
    "datum":"North American Datum 1927"
}'::json)::text ;

result


NOTICE:  DEBUG is on: t
NOTICE:  raw input: {
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"13S",
    "utm_ns":"3885437",
    "utm_ew":"636394",
    "datum":"North American Datum 1927"
}
NOTICE:  converting from UTM SRID
                (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone) 
                to datum SRID (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum)
NOTICE:  utm_srid: 32613
NOTICE:  datum_srid: 4267
NOTICE:  utm_ew: 636394
NOTICE:  utm_ns: 3885437
NOTICE:  conversion operation: ST_Transform(ST_SetSrid(ST_MakePoint(utm_ew,utm_ns),utm_srid),datum_srid)
NOTICE:  ---------------------------------------------------------------------------------
NOTICE:  At this point we should have coordinates transformed to DD.ddd format, and a datum
NOTICE:  noconvlat is not null and noconvlong is not null - we can proceed
NOTICE:  noconvlat: 35.10245630259737
NOTICE:  noconvlong: -103.50299489232316
NOTICE:  v_datum: North American Datum 1927
NOTICE:  datum_srid: 4267
NOTICE:  wgs84: 4326
NOTICE:  datum_srid: 4267
NOTICE:  wgs84: 4326
NOTICE:  Operation: SELECT ST_Transform(ST_SetSrid(ST_MakePoint(noconvlong,noconvlat), datum_srid), wgs84 ) as geom)
                          convertrawcoords                          
--------------------------------------------------------------------
 {"lat":35.102523027164224,"lng":-103.50348741115278,"status":"OK"}
(1 row)
mkoo commented 1 year ago

this all looks good-- however, the SRID for the UTM conversion may only assume values in WGS84 (or similar), so I dont know if this order will work-- maybe transform datum first then convert to DD

let me know when it's hooked up to whatever widget and I can test like a fiend!

I'll say again: I like the idea of a standalone service called Coordinate Converter that we can access whenever

Also, this guy has a nice explanation and recitation on all the math we are trying to avoid with PostGIS (he's the author of the excel calculator): https://stevedutch.net/usefuldata/utmformulas.htm

dustymc commented 1 year ago

I don't know how I'd accomplish "transform datum first then convert to DD" - I have two SRIDs to consider when dealing with UTM, I'm not sure anything other than what I've done could use both?? I'm willing to toss the whole thing if that's what's necessary, I just don't quite see how that could be done - I need more specifics, or some way of understanding.

assume

I have tried to avoid any assumptions at all - if there are any remaining (that don't result in errors) they should be treated as bugs.

standalone service

We've had that for a while, it uses the same code, but it doesn't debug. Edit locality (and probably elsewhere), click the shiny.

Screenshot 2023-03-13 at 4 31 50 PM

If you can get that to return nonsense give me your input and I'll run it through the command-line version with debug on.

explanation

For better or worse, that's all black-box to me - postgis handles that. I have included the operations I've used in the debug, it's definitely possible that I'm making the wrong call, but I suspect spotting that is beyond me.

mkoo commented 1 year ago

Post gis should handle datum transformations out of the box... https://postgis.net/docs/ST_Transform.html

are we barking up the wrong tree? maybe SRID is not the way to go-- PROJ4? (can you check to see if available in our installation: https://postgis.net/docs/PostGIS_Full_Version.html)

On Mon, Mar 13, 2023 at 4:36 PM dustymc @.***> wrote:

I don't know how I'd accomplish "transform datum first then convert to DD"

  • I have two SRIDs to consider when dealing with UTM, I'm not sure anything other than what I've done could use both?? I'm willing to toss the whole thing if that's what's necessary, I just don't quite see how that could be done - I need more specifics, or some way of understanding.

assume

I have tried to avoid any assumptions at all - if there are any remaining (that don't result in errors) they should be treated as bugs.

standalone service

We've had that for a while, it uses the same code, but it doesn't debug. Edit locality (and probably elsewhere), click the shiny.

[image: Screenshot 2023-03-13 at 4 31 50 PM] https://user-images.githubusercontent.com/5720791/224854791-1babd8c3-ae37-4835-ae14-d26e98134db5.png

If you can get that to return nonsense give me your input and I'll run it through the command-line version with debug on.

explanation

For better or worse, that's all black-box to me - postgis handles that. I have included the operations I've used in the debug, it's definitely possible that I'm making the wrong call, but I suspect spotting that is beyond me.

— Reply to this email directly, view it on GitHub https://github.com/ArctosDB/arctos/issues/5910#issuecomment-1467122612, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATH7ULQW6VCLA7IBDF4STDW36VQ5ANCNFSM6AAAAAAVZV575E . You are receiving this because you were mentioned.Message ID: @.***>

dustymc commented 1 year ago

barking up the wrong tree

Could be, I'd not be entirely surprised to find that I'm just not understanding something important and doing it all wrong. (Which is why I've been asking everyone I can think of for help since we got PG!)


arctosprod@arctos>> SELECT PostGIS_Full_Version();
NOTICE:  Topology support cannot be inspected. Is current user granted USAGE on schema "topology" ?
                                                                                 postgis_full_version                                                                                  
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="3.0.3 6660953" [EXTENSION] PGSQL="120" GEOS="3.9.0-CAPI-1.16.2" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" RASTER
(1 row)

Time: 1026.407 ms (00:01.026)

Perhaps pushing TACC on our long-overdue Postgres upgrade makes sense if there's some chance this is related to - well, anything other than my code, I suppose. I've done all I know how to on that front.

mkoo commented 1 year ago

well cant blame it on missing PROJ4 or GDAL at least! It's there

ok let me mull over -- I think basically if the conversion happens in the current order: 1- convert to DD and 2-transform datum; the first will always assume WGS84 and the DD will be in WGS84. and the 2 step is pointless since it's just going from WGS84 to WGS 84

SRID doesnt even have a NAD27/ UTM any zone option. Everything starts and stops at one of the WGS84 ellipsoids

I think the answer lies in how you call up the transformation method/parameters... if you dont use ST_SetSrid what are other options?

dustymc commented 1 year ago

For the UTM (in utm_zone-linked-SRID) conversion:

NOTICE:  conversion operation: ST_Transform(ST_SetSrid(ST_MakePoint(utm_ew,utm_ns),utm_srid),datum_srid)

For the input datum (from ctdatum) to wgs84 conversion, I use

NOTICE:  Operation: SELECT ST_Transform(ST_SetSrid(ST_MakePoint(noconvlong,noconvlat), datum_srid), wgs84 ) as geom)
dustymc commented 1 year ago

IDK if this is easier to read or not:


NOTICE:  conversion operation: 
            ST_Transform(
                ST_SetSrid(
                    ST_MakePoint(
                        utm_ew,
                        utm_ns
                    ),
                    utm_srid
                ),
                datum_srid
            )

and


NOTICE:  Operation: 
            SELECT 
                ST_Transform(
                    ST_SetSrid(
                        ST_MakePoint(
                            noconvlong,
                            noconvlat
                        ),
                        datum_srid
                    ), 
                    wgs84
                )
mkoo commented 1 year ago

I dont know if it's easier, just wavier...

OK here's my current thinking (with encouragement from Deck, confusion all my own) can we try transformation not using SRID but PROJ4 instaed?

See https://postgis.net/docs/ST_Transform.html (and note that the set_SRID is fine for the TO proj but I think given that we cannot find a SRID for NAD 27 then we need to set the geometry in a different way.

right now, the code for the FROM project is NOTICE: utm_srid: 32613 which is https://epsg.io/32613 and is WGS84 we need to start from https://epsg.io/4267

So something like +proj=utm +ellps=clrk66 +nadgrids=NTv2_0.gsb +no_defs +type=crs (probably wrong!)

but look here:https://proj.org/operations/projections/utm.html

BTW, cannot find a postgis UTM conversion tool but plenty of java and even aspx! New Berkeleymapper relies on JTS Topology Suite

dustymc commented 1 year ago

Certainly nothing's off the table, but it seems unlikely to me that our database has forgotten how to define datums?? (And according to the doc that avoids optimization, while I am still struggling with https://github.com/ArctosDB/internal/issues/222.)

mkoo commented 1 year ago

Certainly nothing's off the table, but it seems unlikely to me that our database has forgotten how to define datums??

Not at all what I am suggesting; I'm saying that SRID parameters that are being used do not account for NAD27-- only WGS 84 appears to be available-- someone can help me find it but I cant.

all the references in the SQL to SRID only refers to the WGS84 datum. So maybe try referencing PROJ4 which is installed in our PG (See https://postgis.net/docs/ST_Transform.html ) and has NAD 27 and UTM combos

This is the issue I think:

FIRST OPERATION: convert to DD.dd format using the supplied datum --> the problem is it does not use the supplied datum; the SRID transformation for UTM only assumes WGS84 (see above epsg links)

SECOND OPERATION: at this point we should have DD.dd format data in the supplied datum, just convert it to wgs84 --> needless to say the conversion is pointless since it's going from WGS84 to WGS84

Testing the converter against my GIS seems to bear this out.

dustymc commented 1 year ago

the problem is it does not use the supplied datum;

I don't think that's true. I added a compiled-debug version of the SQL at https://github.com/ArctosDB/PG_DDL/blob/master/function/convertRawCoords.sql#L269

For the example I've been using, it spits out

NOTICE: SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(636394,3885437),32613),4267) as geom) x

copypasta and

arctosprod@arctos>> SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(636394,3885437),32613),4267) as geom) x; st_x | st_y
---------------------+------------------- -103.50299489232316 | 35.10245630259737 (1 row)

636394,3885437 are my UTM coordinates, which I'm making into a point

The srid of that point is 32613 because my input utm zone is https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone#13s

Then I transform that to 4267 because my input datum is https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum#north_american_datum_1927

which, I hope, leaves me with a point in NAD27, which then gets passed on to the WGS84 converter.

I'm pretty clearly lost, but I'm still not seeing any assumptions in there.

mkoo commented 1 year ago

arctosprod@arctos>> SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(636394,3885437),32613),4267) as geom) x;

right there is the assumption: 32613 I believe that is SRID: https://epsg.io/32613 which is "WGS 84 / UTM zone 13N" we need SRID that is NAD 27 / UTM zone 13N which I cannot find (have been complaining see beginning of this!)

UTM conversion has to have the datum stated up front as the math is different if converting from WGS84 UTM and NAD27 UTM. (it only doesnt matter if you are sticking within the ellipsoid family, so WGS84, GRS80, NAD83 essentially all the same)

Maybe I'm not explaining it well-- Make sense? Call tomorrow?

dustymc commented 1 year ago

Yes call, say when. I think maybe I get it - UTM zone and datum are inseparable in PGIS? - but ??

Here's the other example:


select convertRawCoords('{
    "debug":"true",
    "orig_lat_long_units":"decimal degrees",
    "dec_lat":"36.12536",
    "dec_long":"-107.2204214",
    "datum":"North American Datum 1927"
}'::json)::text ;

NOTICE:  DEBUG is on: t
NOTICE:  raw input: {
    "debug":"true",
    "orig_lat_long_units":"decimal degrees",
    "dec_lat":"36.12536",
    "dec_long":"-107.2204214",
    "datum":"North American Datum 1927"
}
NOTICE:  decimal degrees
NOTICE:  ---------------------------------------------------------------------------------
NOTICE:  At this point we should have coordinates transformed to DD.ddd format, and a datum
NOTICE:  noconvlat is not null and noconvlong is not null - we can proceed
NOTICE:  noconvlat: 36.12536
NOTICE:  noconvlong: -107.2204214
NOTICE:  v_datum: North American Datum 1927
NOTICE:  datum_srid: 4267
NOTICE:  wgs84: 4326
NOTICE:  Operation: 
            SELECT 
                ST_Transform(
                    ST_SetSrid(
                        ST_MakePoint(
                            noconvlong,
                            noconvlat
                        ),
                        datum_srid
                    ), 
                    wgs84
                )
                         convertrawcoords                          
-------------------------------------------------------------------
 {"lat":36.12538827808947,"lng":-107.22102918614863,"status":"OK"}
mkoo commented 1 year ago

Holy prime meridians! @dustymc this is what we need to add: list of all SRID values for NAD27 by state plane https://www.spatialreference.org/ref/?page=7&search=nad+27

It took me too long to think of looking stuff up in ESRI files... anyway, starting a list here of what I think we need for UTM conversions

can you test one and then I can update the CT (or you can)

dustymc commented 1 year ago

Oh my.

So https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone is essentially useless, those things can't stand alone.

UTM is completely incompatible with https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum

Correct?

If that's right, then I suggest a not-so-trivial rebuild is necessary

  1. drop utm_zone - it just adds confusion, doesn't supply what we need
  2. flip ctdatum (which is NOT useful as strings, doesn't line up with other stuff, I think results in GBIF "you didn't spell this thing we just make up correctly!!" errors, etc.), so that eg https://epsg.io/4204 (or enough of it to get there) becomes the data, and Ain el Abd 1970 is the metadata. (Optional for my needs, but seems like it would be less-ambiguous for everyone, and would prevent needing to somehow spell out whatever of https://www.spatialreference.org/ref/epsg/26755/ is critical to understanding what "datum" to use)
  3. add whatever we need (eg utm-zone-in-datum) to the one remaining 'srid authority' (whatever we call it) to convert whatever we need to convert.

  1. Am I lost?
  2. Is there a better way?

I don't quite know how to translate our UTM zones to the NAD-having-thingees but I picked one that seemed slightly plausible and this happened:

arctosprod@arctos>>  SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(636394,3885437),26755),4267) as geom) x;
        st_x         |        st_y        
---------------------+--------------------
 -110.91027267851365 | 47.178197285008686

That, I think/hope, is a https://www.spatialreference.org/ref/epsg/26755/ UTM (636394,3885437) reprojected to DD.dd WGS84.

mkoo commented 1 year ago

So https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone is essentially useless, those things can't stand alone.

not useless; it needs to be expanded to include both WGS84 and NAD27 (let's just concentrate on US UTM zones for now)

UTM is completely incompatible with https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum

not incompatible, just not relevant. ctdatum is still needed for converting between datums within DD coordinates

happy to talk this through in realtime ... (seriously this time! just let me know)

mkoo commented 1 year ago

how's this: I can expand the ctutm_zones and correct some of the ones on the list with new values. But you will have to adjust the converter accordingly. Possibly the easier solution is to clean up the dropdown of the utm zones so it also includes the datum-- eg NAD27 / 10N and WGS84 / 10N would both be options. It would make the dropdown longer. Or alternatively you could have two separate fields to chose UTM and then choose zone and leave the the selection of the correct SRID behind the scenes

@dustymc text me if you want to chat

dustymc commented 1 year ago

My concern is that a user can then send

so I think that's functional, and maybe even the best we can do immediately, but it does not seem like a great solution.

@Nicole-Ridgwell-NMMNHS @lin-fred thoughts??

dustymc commented 1 year ago

vague plan

Eventually/maybe: do something unstoopider.

mkoo commented 1 year ago

Please load the attached csv -- that should cover the NAD27 UTM zones for the northern hemisphere ctutm_zone_load.csv

I fixed about half the current table but realized that the MGRS zones really need hemisphere I think to properly determine. Most are not covered by EPSG -- I'll keep looking

Meanwhile once loaded and converter interface updated, let's test the sucker

dustymc commented 1 year ago
ALTER TABLE collecting_event ALTER COLUMN utm_zone TYPE varchar(30);
ALTER TABLE ctutm_zone ALTER COLUMN utm_zone TYPE varchar(30);
ALTER TABLE log_ctutm_zone ALTER COLUMN n_utm_zone TYPE varchar;
ALTER TABLE log_ctutm_zone ALTER COLUMN o_utm_zone TYPE varchar;
insert into ctutm_zone(utm_zone,description,srid)(select utm_zone,description,srid::int  from temp_cache.temp_dlm_uptbl );
update ctutm_zone set utm_zone =trim(utm_zone);

And the UTM from above fails as it should:


select convertRawCoords('{
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"13S",
    "utm_ns":"3885437",
    "utm_ew":"636394",
    "datum":"North American Datum 1927"
}'::json)::text ;

NOTICE:  raw input: {
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"13S",
    "utm_ns":"3885437",
    "utm_ew":"636394",
    "datum":"North American Datum 1927"
}
NOTICE:  converting from UTM SRID
                (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone) 
                to datum SRID (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum)
                            convertrawcoords                            
------------------------------------------------------------------------
 {"status":"fail","message":"incompatible datum","lat":null,"lng":null}
(1 row)

and with a new zone:


NOTICE:  DEBUG is on: t
NOTICE:  raw input: {
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"NAD27 / UTM zone 13N",
    "utm_ns":"3885437",
    "utm_ew":"636394",
    "datum":"North American Datum 1927"
}
NOTICE:  converting from UTM SRID
                (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone) 
                to datum SRID (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum)
NOTICE:  utm_srid: 26713
NOTICE:  datum_srid: 4267
NOTICE:  utm_ew: 636394
NOTICE:  utm_ns: 3885437
NOTICE:  conversion operation: 
            ST_Transform(
                ST_SetSrid(
                    ST_MakePoint(
                        utm_ew,
                        utm_ns
                    ),
                    utm_srid
                ),
                datum_srid
            )
NOTICE:   SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(636394,3885437),26713),4267) as geom) x
NOTICE:  ---------------------------------------------------------------------------------
NOTICE:  At this point we should have coordinates transformed to DD.ddd format, and a datum
NOTICE:  noconvlat is not null and noconvlong is not null - we can proceed
NOTICE:  noconvlat: 35.10430694762628
NOTICE:  noconvlong: -103.50348950441393
NOTICE:  v_datum: North American Datum 1927
NOTICE:  datum_srid: 4267
NOTICE:  wgs84: 4326
NOTICE:  Operation: 
            SELECT 
                ST_Transform(
                    ST_SetSrid(
                        ST_MakePoint(
                            noconvlong,
                            noconvlat
                        ),
                        datum_srid
                    ), 
                    wgs84
                )
                          convertrawcoords                          
--------------------------------------------------------------------
 {"lat":35.104373629973686,"lng":-103.50398204883706,"status":"OK"}

If that's right let's CSV all new datums - I've got the dependency hardcoded in (https://github.com/ArctosDB/PG_DDL/blob/master/function/convertRawCoords.sql#L252), the zone needs to consist of NAD27 or WGS84 and then a slash, accompanied by the right datum, or it'll fail.

mkoo commented 1 year ago

on test, right? let me try it but seems like that works!

dustymc commented 1 year ago

@Nicole-Ridgwell-NMMNHS @lin-fred @mkoo can I get some examples to test this with? If this is solved I need to get a CSV of the new zones back to production, if it's not then I'd like to know ASAP.

I tried to use https://github.com/ArctosDB/arctos/issues/5279 but I don't know what zone it should be - "13N" doesn't seem to be a viable possibility??

Nicole-Ridgwell-NMMNHS commented 1 year ago

Thanks @mkoo @dustymc for working on this. I'm catching up from being out last week, but I'll see if I can put together a set of examples by the end of the day. #5279 should work, it is northern hemisphere, you can use 13N.

dustymc commented 1 year ago

Thanks @Nicole-Ridgwell-NMMNHS . We need documentation - at some point I'd convinced myself that "S" is a zone, not a direction, but now a couple comments make me think that's not right, or too simple, or ???

https://www.dmap.co.uk/utmworld.htm

Screenshot 2023-03-20 at 9 59 59 AM

If that's right then the below is wrong - I used 13N (Colombia-ish) instead of 13S (New Mexico-ish, and not yet in the code table).

Anyway...here's a type-from-image attempt at 5279, is this doing what it should be (assuming I didn't fatfinger anything!)?


select convertRawCoords('{
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"NAD27 / UTM zone 13N",
    "utm_ns":"3988678",
    "utm_ew":"285729",
    "datum":"North American Datum 1927"
}'::json)::text ;

NOTICE:  DEBUG is on: t
NOTICE:  raw input: {
    "debug":"true",
    "orig_lat_long_units":"UTM",
    "utm_zone":"NAD27 / UTM zone 13N",
    "utm_ns":"3988678",
    "utm_ew":"285729",
    "datum":"North American Datum 1927"
}
NOTICE:  converting from UTM SRID
                (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctutm_zone) 
                to datum SRID (https://arctos.database.museum/info/ctDocumentation.cfm?table=ctdatum)
NOTICE:  utm_srid: 26713
NOTICE:  datum_srid: 4267
NOTICE:  utm_ew: 285729
NOTICE:  utm_ns: 3988678
NOTICE:  conversion operation: 
            ST_Transform(
                ST_SetSrid(
                    ST_MakePoint(
                        utm_ew,
                        utm_ns
                    ),
                    utm_srid
                ),
                datum_srid
            )
NOTICE:   SELECT ST_X(geom), ST_Y(geom) from (SELECT ST_Transform(ST_SetSrid(ST_MakePoint(285729,3988678),26713),4267) as geom) x
NOTICE:  ---------------------------------------------------------------------------------
NOTICE:  At this point we should have coordinates transformed to DD.ddd format, and a datum
NOTICE:  noconvlat is not null and noconvlong is not null - we can proceed
NOTICE:  noconvlat: 36.02086763136399
NOTICE:  noconvlong: -107.37777549614447
NOTICE:  v_datum: North American Datum 1927
NOTICE:  datum_srid: 4267
NOTICE:  wgs84: 4326
NOTICE:  Operation: 
            SELECT 
                ST_Transform(
                    ST_SetSrid(
                        ST_MakePoint(
                            noconvlong,
                            noconvlat
                        ),
                        datum_srid
                    ), 
                    wgs84
                )
                          convertrawcoords                          
--------------------------------------------------------------------
 {"lat":36.020897237216225,"lng":-107.37838703162896,"status":"OK"}
Nicole-Ridgwell-NMMNHS commented 1 year ago

The latitude bands (C-M southern hemisphere and N-X northern hemisphere) are kind of an addition to UTM (see https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Latitude_bands) and aren't necessary for conversion. For conversion, all that matters is northern or southern hemisphere. For example, see the NCAT converter https://www.ngs.noaa.gov/NCAT/ (select UTM).

Nicole-Ridgwell-NMMNHS commented 1 year ago

Here is a file of UTMs to test. Incorrect_DEC_LAT/LONG are what bulk edit locality produced from the UTMs when I loaded the UTM data last year. Correct_DEC_LAT/LONG are UTM conversions from the NCAT converter. All zones are northern hemisphere.

TestUTMs.xlsx

dustymc commented 1 year ago

I did this:

create table temp_test_cc as select * from temp_cache.temp_dlm_uptbl ;
alter table temp_test_cc add dlat numeric;
alter table temp_test_cc add dlon numeric;
alter table temp_test_cc add utmzone varchar;
alter table temp_test_cc add key serial not null;

update temp_test_cc set utmzone='NAD27 / UTM zone 13N' where zone='13';
update temp_test_cc set utmzone='NAD27 / UTM zone 18N' where zone='18';
update temp_test_cc set utmzone='NAD27 / UTM zone 12N' where zone='12';
update temp_test_cc set utmzone='NAD27 / UTM zone 16N' where zone='16';
update temp_test_cc set utmzone='NAD27 / UTM zone 15N' where zone='15';
update temp_test_cc set utmzone='NAD27 / UTM zone 14N' where zone='14';

CREATE OR REPLACE function test_test() returns void AS $body$
DECLARE
c int;
r record;
la numeric;
lo numeric;
jr json;
BEGIN
for r in (select * from temp_test_cc) loop
select convertRawCoords(json_build_object(
      'debug','false',
    'orig_lat_long_units','UTM',
    'datum',r.datum,
    'utm_zone',r.utmzone,
    'utm_ew',r.easting,
    'utm_ns',r.northing
    )) into jr;

la=jr->>'lat';
lo=jr->>'lng';
update temp_test_cc set dlat=la,dlon=lo where key=r.key;
end loop;
END;
$body$
LANGUAGE PLPGSQL
SECURITY DEFINER
;
select test_test();

result:

temp_test_cc.csv

yay??

Nicole-Ridgwell-NMMNHS commented 1 year ago

Looks good to me!

dustymc commented 1 year ago

@mkoo do you have csv I can use for the code table in prod, or should I try to patch that together, or ???

dustymc commented 1 year ago

I brought the code table data from test over to prod, and then I cleaned up and inserted the documentation as new values. There's currently some redundancy in the table, but it's not functional redundancy (the old values will throw errors even though they share SRID with values that will convert) and I couldn't see any other way to make sure the maybe-buggy-depending-on-path data don't get mixed up with the zone-and-datum converted data.

Tentatively closing this, I think we're done here.

mkoo commented 1 year ago

sorry so late for all these issues-- did you get my csv from the other issue? It looks like the other NAD27 flavors did not get loaded. Let me file a new issue and send along!

dustymc commented 1 year ago

the UI caches, here's current prod data:

arctosprod@arctos>> 
arctosprod@arctos>>  select * from ctutm_zone order by utm_zone;
       utm_zone       |      description      | srid  
----------------------+-----------------------+-------
 10N                  | WGS 84 / UTM zone 10N | 32610
 10S                  | WGS 84 / UTM zone 10S | 32710
 10T                  | 10T                   | 32610
 11N                  | WGS 84 / UTM zone 11N | 32611
 11S                  | WGS 84 / UTM zone 11S | 32711
 11T                  | 11T (MGRS)            | 32611
 12N                  | WGS 84 / UTM zone 12N | 32612
 12R                  | 12R (MGRS)            | 32612
 12S                  | WGS 84 / UTM zone 12S | 32712
 12T                  | 12T (MGRS)            | 32612
 13E                  | 13E (MGRS)            | 32713
 13N                  | WGS 84 / UTM zone 13N | 32613
 13R                  | 13R (MGRS)            | 32613
 13S                  | WGS 84 / UTM zone 13S | 32713
 13T                  | 13T (MGRS)            | 32613
 14N                  | WGS 84 / UTM zone 14N | 32614
 14R                  | 14R (MGRS)            | 32614
 14S                  | WGS 84 / UTM zone 14S | 32714
 15N                  | WGS 84 / UTM zone 15N | 32615
 15S                  | WGS 84 / UTM zone 15S | 32715
 15T                  | 15T (MGRS)            | 32615
 15U                  | 15U (MGRS)            | 32615
 16N                  | WGS 84 / UTM zone 16N | 32616
 16S                  | WGS 84 / UTM zone 16S | 32716
 16T                  | 16T                   | 32616
 16U                  | 16U                   | 32616
 17P                  | 17P                   | 32617
 17S                  | 17S                   | 32617
 17T                  | 17T                   | 32617
 17U                  | 17U                   | 32617
 18T                  | 18T                   | 32618
 19S                  | 19S                   | 32619
 36K                  | 36K                   | 32736
 36L                  | 36L                   | 32736
 37K                  | 37K                   | 32737
 37L                  | 37L                   | 32737
 48P                  | 48P                   | 32648
 5W                   | 5W                    |  3751
 6N                   | WGS 84 / UTM zone 6N  |  3713
 6V                   | 6V                    |  3713
 6W                   | 6W                    |  3713
 NAD27 / UTM zone 10N |  NAD27 / UTM zone 10N | 26710
 NAD27 / UTM zone 11N |  NAD27 / UTM zone 11N | 26711
 NAD27 / UTM zone 12N |  NAD27 / UTM zone 12N | 26712
 NAD27 / UTM zone 13N |  NAD27 / UTM zone 13N | 26713
 NAD27 / UTM zone 14N |  NAD27 / UTM zone 14N | 26714
 NAD27 / UTM zone 15N |  NAD27 / UTM zone 15N | 26715
 NAD27 / UTM zone 16N |  NAD27 / UTM zone 16N | 26716
 NAD27 / UTM zone 17N |  NAD27 / UTM zone 17N | 26717
 NAD27 / UTM zone 18N |  NAD27 / UTM zone 18N | 26718
 NAD27 / UTM zone 19N |  NAD27 / UTM zone 19N | 26719
 NAD27 / UTM zone 1N  | NAD27 / UTM zone 1N   | 26701
 NAD27 / UTM zone 20N |  NAD27 / UTM zone 20N | 26720
 NAD27 / UTM zone 21N |  NAD27 / UTM zone 21N | 26721
 NAD27 / UTM zone 22N |  NAD27 / UTM zone 22N | 26722
 NAD27 / UTM zone 2N  |  NAD27 / UTM zone 2N  | 26702
 NAD27 / UTM zone 3N  |  NAD27 / UTM zone 3N  | 26703
 NAD27 / UTM zone 4N  |  NAD27 / UTM zone 4N  | 26704
 NAD27 / UTM zone 5N  |  NAD27 / UTM zone 5N  | 26705
 NAD27 / UTM zone 6N  |  NAD27 / UTM zone 6N  | 26706
 NAD27 / UTM zone 7N  |  NAD27 / UTM zone 7N  | 26707
 NAD27 / UTM zone 8N  |  NAD27 / UTM zone 8N  | 26708
 NAD27 / UTM zone 9N  |  NAD27 / UTM zone 9N  | 26709
 WGS84 / UTM zone 10N | WGS84 / UTM zone 10N  | 32610
 WGS84 / UTM zone 10S | WGS84 / UTM zone 10S  | 32610
 WGS84 / UTM zone 11N | WGS84 / UTM zone 11N  | 32611
 WGS84 / UTM zone 11S | WGS84 / UTM zone 11S  | 32711
 WGS84 / UTM zone 12N | WGS84 / UTM zone 12N  | 32612
 WGS84 / UTM zone 12S | WGS84 / UTM zone 12S  | 32712
 WGS84 / UTM zone 13N | WGS84 / UTM zone 13N  | 32613
 WGS84 / UTM zone 13S | WGS84 / UTM zone 13S  | 32713
 WGS84 / UTM zone 14N | WGS84 / UTM zone 14N  | 32614
 WGS84 / UTM zone 14S | WGS84 / UTM zone 14S  | 32714
 WGS84 / UTM zone 15N | WGS84 / UTM zone 15N  | 32615
 WGS84 / UTM zone 15S | WGS84 / UTM zone 15S  | 32715
 WGS84 / UTM zone 16N | WGS84 / UTM zone 16N  | 32616
 WGS84 / UTM zone 16S | WGS84 / UTM zone 16S  | 32716
 WGS84 / UTM zone 6N  | WGS84 / UTM zone 6N   |  3713
(78 rows)
mkoo commented 1 year ago

oh shoot! Ignore #5279 !

You can see that we still have a little more clean up on the other subzones like K-T's -- can you let us know the datum needed? for now, you can use the closest 'parent' utm zone. (sorry that's not their actual names but getting low on blood sugar here....)

dustymc commented 1 year ago

Not sure I'm following, ditto on the blood sugar, but here's my code:


        if not (v_datum='North American Datum 1927' and trim(split_part(v_utm_zone,'/',1))='NAD27') or 
            (v_datum ='World Geodetic System 1984' and trim(split_part(v_utm_zone,'/',1))='WGS84') then
            select 
                'fail' as status,

so right now I can only deal with SRIDs that correspond to WGS84 or NAD27, and the utm_zone values needs to start with those then space-slash then whatever. More zones in those datums is no problem, just be consistent in the spelling, and for other datums I'll need a similar pattern and a quick rewrite.

Nicole-Ridgwell-NMMNHS commented 1 year ago

Coordinate converter doesn't seem to be working for WGS 84 data:

image

dustymc commented 1 year ago

Try now

Nicole-Ridgwell-NMMNHS commented 1 year ago

Thanks, that fixed it!