Merge pull request #112 from jlaura/proj

# Changelog

All changes that impact users of this project will be documented in this file.

The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),

and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

<!---

This document is intended for users of the applications and API. Changes to things

like tests should not be noted in this document.

When updating this file for a PR, add an entry for your change under Unreleased

and one of the following headings:

 - Added - for new features.

 - Changed - for changes in existing functionality.

 - Deprecated - for soon-to-be removed features.

 - Removed - for now removed features.

 - Fixed - for any bug fixes.

 - Security - in case of vulnerabilities.

If the heading does not yet exist under Unreleased, then add it as a 3rd heading,

with three #.

When preparing for a public release candidate add a new 2nd heading, with two #, under

Unreleased with the version number and the release date, in year-month-day

format. Then, add a link for the new version at the bottom of this document and

update the Unreleased link so that it compares against the latest release tag.

When preparing for a bug fix release create a new 2nd heading above the Fixed

heading to indicate that only the bug fixes and security fixes are in the bug fix

release.

-->

## Unreleased

### Changed

- Removed all `pyproj` calls from csm.py, abstracting them into the reprojection and pyproj.Transformer code inside utils.py. Updated the transformations to use the new pipeline style syntax to avoid deprecation warnings about old syntax.p

import jinja2

from osgeo import ogr

import numpy as np

import pyproj

import requests

import scipy.stats

from functools import singledispatch

from plio.io.io_gdal import GeoDataset

from knoten import utils

class NumpyEncoder(json.JSONEncoder):

    def default(self, obj):

        if isinstance(obj, np.ndarray):

    intersection = generate_ground_point(0.0, image_pt, camera)

    iterations = 0

    semi_major, semi_minor = get_radii(camera)

    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)

    source_proj = f'+proj=cart +a={semi_major} +b={semi_minor}'

    dest_proj = f'+proj=lonlat +a={semi_major} +b={semi_minor}'

    transformer = utils.create_transformer(source_proj, dest_proj)

    while iterations != max_its:

        lon, lat, alt = pyproj.transform(ecef, lla, intersection.x, intersection.y, intersection.z)

        lon, lat, alt = transformer.transform(intersection.x, 

                                              intersection.y, 

                                              intersection.z,

                                              errcheck=True)

        px, py = dem.latlon_to_pixel(lat, lon)

        height = dem.read_array(1, [px, py, 1, 1])[0][0]

    else:

        semi_major, semi_minor = radii

    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)

    source_proj = f'+proj=cart +a={semi_major} +b={semi_minor}'

    dest_proj = f'+proj=lonlat +a={semi_major} +b={semi_minor}'

    transformer = utils.create_transformer(source_proj, dest_proj)

    gnds = np.empty((len(boundary), 3))

        gnds[i] = [gnd.x, gnd.y, gnd.z]

    lons, lats, alts = pyproj.transform(ecef, lla, gnds[:,0], gnds[:,1], gnds[:,2])

    lons, lats, alts = transformer.transform(gnds[:,0], gnds[:,1], gnds[:,2])

    return lons, lats, alts

def generate_gcps(camera, boundary, radii=None):

    latlon_fp = generate_latlon_footprint(camera, boundary, radii=radii)

    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)

    source_proj = f'+proj=lonlat +a={semi_major} +b={semi_minor}'

    dest_proj = f'+proj=cart +a={semi_major} +b={semi_minor}'

    transformer = utils.create_transformer(source_proj, dest_proj)

    # Step over all geometry objects in the latlon footprint

    for i in range(latlon_fp.GetGeometryCount()):

        latlon_coords = np.array(latlon_fp.GetGeometryRef(i).GetGeometryRef(0).GetPoints())

        # Check if the geometry object is populated with points

        if len(latlon_coords) > 0:

            x, y, z = pyproj.transform(lla, ecef,  latlon_coords[:,0], latlon_coords[:,1], latlon_coords[:,2])

            x, y, z = transformer.transform(latlon_coords[:,0], latlon_coords[:,1], latlon_coords[:,2])

            # Step over all coordinate points in a geometry object and update said point

            for j, _ in enumerate(latlon_coords):

    Transformed coordinates as y, x, z

    """

    source_pyproj = pyproj.Proj(proj = source_proj, a = semi_major, b = semi_minor)

    dest_pyproj = pyproj.Proj(proj = dest_proj, a = semi_major, b = semi_minor)

    y, x, z = pyproj.transform(source_pyproj, dest_pyproj, record[0], record[1], record[2], **kwargs)

    transformer = pyproj.Transformer.from_crs(f'+proj={source_proj} +a={semi_major} +b={semi_minor}',

                                              f'+proj={dest_proj} +a={semi_major} +b={semi_minor}',

                                              always_xy=True)

    source_proj = f'+proj={source_proj} +a={semi_major} +b={semi_minor}'

    dest_proj = f'+proj={dest_proj} +a={semi_major} +b={semi_minor}'

    transformer = create_transformer(source_proj, dest_proj)

    x, y, z = transformer.transform(record[0], record[1], record[2], errcheck=True)

    return y, x, z

def create_transformer(source_proj, dest_proj):

    return pyproj.Transformer.from_crs(source_proj,

                                       dest_proj,

                                       always_xy=True)

 No newline at end of file

from knoten import utils

def test_reproject():

    with mock.patch('pyproj.transform', return_value=[1,1,1]) as mock_pyproj:

        res = utils.reproject([1,1,1], 10, 10, 'geocent', 'latlon')

        mock_pyproj.assert_called_once()

        assert res == (1,1,1)

    res = utils.reproject([0,1,0], 10, 10, 'geocent', 'latlon')

    print(res)

    assert res[0] == 0

    assert res[1] == 90

    assert res[2] == -9