Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

diff/show: Add reprojection #213

Merged
merged 3 commits into from
Aug 20, 2020
Merged

diff/show: Add reprojection #213

merged 3 commits into from
Aug 20, 2020

Conversation

craigds
Copy link
Member

@craigds craigds commented Aug 18, 2020

Description

diff/show: Add reprojection using the --crs option.

The form of the argument is quite flexible. Equivalent forms:

  • --crs 2193

  • --crs epsg:2193

  • --crs EPSG:2193

  • --crs '+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'

  • --crs 'PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","2193"]]'

  • --crs @./nztm.wkt

  • --crs @./nztm.prj

I have not added any handling for cases where reprojection fails. The process will error (RuntimeError normally) in that case.

Related links:

none

Checklist:

  • Have you reviewed your own change?
  • Have you included test(s)?
  • Have you updated the changelog?

@craigds craigds requested review from olsen232 and rcoup August 18, 2020 00:32
tests/test_diff.py Outdated Show resolved Hide resolved
sno/diff.py Outdated Show resolved Hide resolved
@rcoup
Copy link
Member

rcoup commented Aug 18, 2020

Can we use crs instead of srs (everywhere)? Seems like OGC is standardising on CRS now, no reason we shouldn't as well.

@rcoup
Copy link
Member

rcoup commented Aug 18, 2020

So, these are the formats proj_create() accepts, I don't think we need to support them all right now, but I don't think we should add anything not in this list:

  1. a proj-string,
  2. a WKT string,
  3. an object code (like “EPSG:4326”, “urn:ogc:def:crs:EPSG::4326”, “urn:ogc:def:coordinateOperation:EPSG::1671”),
  4. an Object name. e.g “WGS 84”, “WGS 84 / UTM zone 31N”. In that case as uniqueness is not guaranteed, heuristics are applied to determine the appropriate best match.
  5. a OGC URN combining references for compound coordinate reference systems (e.g “urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717” or custom abbreviated syntax “EPSG:2393+5717”),
  6. a OGC URN combining references for references for projected or derived CRSs e.g. for Projected 3D CRS “UTM zone 31N / WGS 84 (3D)”: “urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031” (added in 6.2)
  7. a OGC URN combining references for concatenated operations (e.g. “urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618”)
  8. a PROJJSON string. The jsonschema is at https://proj.org/schemas/v0.2/projjson.schema.json (added in 6.2)
  9. a compound CRS made from two object names separated with ” + “. e.g. “WGS 84 + EGM96 height” (added in 7.1)

  • --srs 2193

I'm not keen on this one, c/- the above logic.

  • --srs EPSG:2193
  • --srs EPSG:2193

Meant to be the same?

  • --srs '+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'

Presumably this is any proj-string, so eg: --srs '+init=epsg:2193' works too?

Idle not-for-now query: if I add a +nadgrids=/path/to/some.gsb presumably that isn't helpful to anyone. I guess this goes away now new proj basically has every grid in it and prefers them.

--srs @./nztm.wkt
--srs @./nztm.prj

Can we make the content of this just parsed like the command-line argument? Ie could contain a proj string, EPSG:1234, or anything else we support?

And related to axis-order handling ideas, can we add URN support? --crs=urn:ogc:def:crs:EPSG::4326

Copy link
Member

@rcoup rcoup left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding some very basic handling (ie. error message + quit) from reprojection failures would be good.

sno/diff.py Outdated Show resolved Hide resolved
sno/diff.py Outdated Show resolved Hide resolved
sno/diff.py Outdated Show resolved Hide resolved
sno/diff.py Outdated Show resolved Hide resolved
@@ -392,19 +408,23 @@ def _out(dataset, diff):
class LazyGeojsonFeatureOutput:
"""Wrapper of KeyValue that lazily serialises it as GEOJSON when sent to json.dumps"""

__slots__ = ("change_type", "key_value")
__slots__ = ("change_type", "key_value", "geometry_transform")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

isn't there one transform per dataset? Do we need to carry it around with every feature?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, there's only one of these objects (all features get the same one).

but the geometry isn't turned into an OGR geometry until it's being streamed out as JSON. So we need to pass the transform down into the lazily-serialized object so it can be used after we have an OGR geometry.

Doing it this way is quite elegant because it ensures minimal memory overhead; we don't have to store millions of OGR geometries in memory (we just create one at a time while streaming JSON out)

Copy link
Member

@rcoup rcoup Aug 19, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sure, doing the transform at emit-time makes total sense. But why do we need to pass N million pointers to the same Transform object about with each feature?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the alternative being a global variable? I try to avoid those where it's not absolutely necessary

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Create an output context containing that?

I'm not familiar enough with how this is working though, I'd have to dig in tomorrow — eg: I thought we were generating/producing a feature at a time, but that doesn't seem to be the case, we're building up big features lists.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sort of - we're building up lists of dicts, but only transforming & converting during json.dump(). It's a lot less work up-front than it used to be.

sno/diff_output.py Show resolved Hide resolved
sno/geometry.py Outdated Show resolved Hide resolved
sno/diff_output.py Outdated Show resolved Hide resolved
@craigds
Copy link
Member Author

craigds commented Aug 18, 2020

  • --srs 2193

I'm not keen on this one, c/- the above logic.

ok, will remove

Presumably this is any proj-string, so eg: --srs '+init=epsg:2193' works too?

yes it does

Idle not-for-now query: if I add a +nadgrids=/path/to/some.gsb presumably that isn't helpful to anyone. I guess this goes away now new proj basically has every grid in it and prefers them.

yeah, that'll cause issues, but let's not deal with it. User problem

--srs @./nztm.wkt
--srs @./nztm.prj
Can we make the content of this just parsed like the command-line argument? Ie could contain a proj string, EPSG:1234, or anything else we support?

it is - I wasn't suggested the filename here is important. It just reads the file and parses it as a string like everything else

And related to axis-order handling ideas, can we add URN support? --crs=urn:ogc:def:crs:EPSG::4326

will see if I can get GDAL to handle them. Otherwise I'm not terribly familiar with them and don't really want to figure out how to parse them correctly

@craigds craigds force-pushed the diff-reprojection branch 2 times, most recently from 1bcbd81 to 7d9cdff Compare August 19, 2020 02:05
@rcoup
Copy link
Member

rcoup commented Aug 19, 2020

And related to axis-order handling ideas, can we add URN support? --crs=urn:ogc:def:crs:EPSG::4326

will see if I can get GDAL to handle them. Otherwise I'm not terribly familiar with them and don't really want to figure out how to parse them correctly

SetFromUserInput() handles them, so done?

@craigds craigds force-pushed the diff-reprojection branch 2 times, most recently from cae6bbf to 53a7447 Compare August 20, 2020 21:09
@craigds craigds requested a review from rcoup August 20, 2020 21:39
@craigds craigds merged commit 43168f8 into master Aug 20, 2020
@craigds craigds deleted the diff-reprojection branch August 20, 2020 23:44
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants