geowatch.gis.spatial_reference module¶

Tools relating to coordinate reference systems

Notes

GEOJSON is assumed to be specified as reversed WGS84 with reversed axes (i.e. lon/lat) unless another CRS is specified.
WGS84 (aka EPSG 4326) is always given as latitude longitude
EPSG stands for “European Petroleum Survey Group”

class geowatch.gis.spatial_reference.RPCTransform(rpcs, elevation='gtop30', axis_mapping='OAMS_TRADITIONAL_GIS_ORDER')[source]¶

Bases: object

Wrapper around rasterio RPC transforms

Parameters:

rpcs (rasterio.rpc.RPC) – rasterio RPC data
elevation (str) – method used to determine the elevation when RPC information is used. Available options are:
- “open-elevation”
- “gtop30”

Notes

By definition rpcs are always referenced against WGS84 (EPSG:4326) coordinates.
However, we are accept and return lon/lat (i.e. reversed WGS84) points here.
TODO: don’t use reversed. Use authority compliant

References

https://rasterio.readthedocs.io/en/latest/topics/reproject.html#reprojecting-with-other-georeferencing-metadata http://geotiff.maptools.org/rpc_prop.html https://github.com/mapbox/rasterio/blob/master/rasterio/rpc.py https://github.com/mapbox/rasterio/blob/master/rasterio/_transform.pyx https://gdal.org/doxygen/gdal__alg_8h.html#af4c3c0d4c79218995b3a1f0bac3700a0

Example

>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> from geowatch.gis.spatial_reference import *  # NOQA
>>> self = RPCTransform.demo()
>>> pxl_pts = np.array([
>>>     [    0,     0],
>>>     [    0, 20000],
>>>     [20000,     0],
>>>     [20000, 20000],
>>> ])
>>> wld_pts = self.warp_world_from_pixel(pxl_pts)
>>> print('wld_pts =\n{}'.format(ub.urepr(wld_pts, nl=1, precision=2)))
wld_pts =
np.array([[128.87,  37.8 ],
          [128.87,  37.7 ],
          [129.  ,  37.81],
          [129.  ,  37.71]]...)

>>> wld_pts = np.array([[self.rpcs.long_off, self.rpcs.lat_off]])
>>> pxl_pts = self.warp_pixel_from_world(wld_pts)
>>> print('pxl_pts = {}'.format(ub.urepr(pxl_pts, nl=1, precision=2)))
pxl_pts = np.array([[17576.2 , 10690.73]]...)

classmethod demo(**kw)[source]¶: Return a demo RPCTransform for testing purposes

classmethod from_gdal(rpc_info, **kwargs)[source]¶

Build an RPC transform from GDAL-style metadata

Parameters:

rpc_info (dict) – gdal RPC dictionary. Typically aquired from gdal.Open(<path>).GetMetadata(domain='RPC')
**kwargs – passed to class:RPCTransform

warp_pixel_from_world(pts_in, return_elevation=False)[source]¶

Parameters:: pts_in (ndarray) – Either an Nx2 array of lon/lat WGS84 coordinates or an Nx3 array of lon/lat/elevation coordinates.

Todo

[ ] Handle CRS

Returns:: An Nx2 array of pixel coordinates
Return type:: pts_out (ndarray)

warp_world_from_pixel(pts_in, return_elevation=False)[source]¶

Parameters:: pts_in (ndarray) – Either an Nx2 array of x,y pixel coordinates, or an Nx3 array of x,y,elevation coordinates.

Todo

[ ] Handle CRS

Returns:: An Nx2 array of lon/lat WGS84 coordinates
Return type:: pts_out (ndarray)

Example

>>> # xdoctest: +SKIP
>>> # xdoctest: +REQUIRES(env:SLOW_DOCTEST)
>>> from geowatch.gis.spatial_reference import *  # NOQA
>>> self = RPCTransform.demo(elevation='gtop30')
>>> pts_in = pxl_pts = np.array([
>>>     [    0,     0],
>>>     [    0, 20000],
>>>     [20000,     0],
>>>     [20000, 20000],
>>> ])
>>> wld_pts = self.warp_world_from_pixel(pxl_pts)
>>> print('wld_pts =\n{}'.format(ub.urepr(wld_pts, nl=1, precision=2)))
>>> #
>>> import osgeo
>>> from geowatch.gis.spatial_reference import *  # NOQA
>>> gpath = '/home/joncrall/data/dvc-repos/smart_watch_dvc/drop0/KR-Pyeongchang-WV/_assets/20140131_a_KRG_011778204_10_0/011778204010_01_003/011778204010_01/011778204010_01_P002_PAN/14JAN31020440-P1BS-011778204010_01_P002.NTF'
>>> #gpath = '/home/joncrall/data/dvc-repos/smart_watch_dvc/drop0/KR-Pyeongchang-WV/_assets/20140131_a_KRG_011778204_10_0/011778204010_01_003/011778204010_01/011778204010_01_P001_PAN/14JAN31020439-P1BS-011778204010_01_P001.NTF'
>>> from osgeo import gdal
>>> ref = gdal.Open(gpath)
>>> rpc_info = ref.GetMetadata(domain='RPC')
>>> pts_in = pxl_pts = np.array([
>>>     [    0,                         0],
>>>     [    0,           ref.RasterYSize],
>>>     [ref.RasterXSize, ref.RasterYSize],
>>>     [ref.RasterXSize,               0],
>>> ])
>>> self = RPCTransform.from_gdal(rpc_info, elevation='gtop30')
>>> wld_pts = self.warp_world_from_pixel(pxl_pts)
>>> kwimage.Polygon(exterior=wld_pts).draw(color='purple', alpha=0.5)
>>> #
>>> self = RPCTransform.from_gdal(rpc_info, elevation='open-elevation')
>>> self.warp_world_from_pixel(pxl_pts)

make_warp_pixel_from_world()[source]¶: Hack for pickelability

make_warp_world_from_pixel()[source]¶: Hack for pickelability

class geowatch.gis.spatial_reference.RPCPixelFromWorldTransform(rpc_transform)[source]¶

Bases: object

Helper class to pickle the warp_pixel_from_world method. I’m not sure if this is really needed. I would think a method can be pickled if its class can be pickeld…

class geowatch.gis.spatial_reference.RPCWorldFromPixelTransform(rpc_transform)[source]¶

Bases: object

Helper class to pickle the warp_world_from_pixel method. I’m not sure if this is really needed. I would think a method can be pickled if its class can be pickeld…