include "base.pxi" cimport cython from cpython.mem cimport PyMem_Free, PyMem_Malloc import copy import re import warnings from collections import namedtuple from pyproj._compat cimport cstrencode from pyproj._context cimport _clear_proj_error, _get_proj_error, pyproj_context_create from pyproj._crs cimport ( _CRS, Base, CoordinateOperation, _get_concatenated_operations, _to_proj4, _to_wkt, create_area_of_use, ) from pyproj._context import _LOGGER, get_context_manager from pyproj.aoi import AreaOfInterest from pyproj.enums import ProjVersion, TransformDirection, WktVersion from pyproj.exceptions import ProjError _AUTH_CODE_RE = re.compile(r"(?P\w+)\:(?P\w+)") cdef dict _TRANSFORMER_TYPE_MAP = { PJ_TYPE_UNKNOWN: "Unknown Transformer", PJ_TYPE_CONVERSION: "Conversion Transformer", PJ_TYPE_TRANSFORMATION: "Transformation Transformer", PJ_TYPE_CONCATENATED_OPERATION: "Concatenated Operation Transformer", PJ_TYPE_OTHER_COORDINATE_OPERATION: "Other Coordinate Operation Transformer", } Factors = namedtuple( "Factors", [ "meridional_scale", "parallel_scale", "areal_scale", "angular_distortion", "meridian_parallel_angle", "meridian_convergence", "tissot_semimajor", "tissot_semiminor", "dx_dlam", "dx_dphi", "dy_dlam", "dy_dphi", ], ) Factors.__doc__ = """ .. versionadded:: 2.6.0 These are the scaling and angular distortion factors. See PROJ :c:type:`PJ_FACTORS` documentation. Parameters ---------- meridional_scale: list[float] Meridional scale at coordinate. parallel_scale: list[float] Parallel scale at coordinate. areal_scale: list[float] Areal scale factor at coordinate. angular_distortion: list[float] Angular distortion at coordinate. meridian_parallel_angle: list[float] Meridian/parallel angle at coordinate. meridian_convergence: list[float] Meridian convergence at coordinate. Sometimes also described as *grid declination*. tissot_semimajor: list[float] Maximum scale factor. tissot_semiminor: list[float] Minimum scale factor. dx_dlam: list[float] Partial derivative of coordinate. dx_dphi: list[float] Partial derivative of coordinate. dy_dlam: list[float] Partial derivative of coordinate. dy_dphi: list[float] Partial derivative of coordinate. """ cdef PJ_DIRECTION get_pj_direction(object direction) except *: # optimized lookup to avoid creating a new instance every time # gh-1205 if not isinstance(direction, TransformDirection): direction = TransformDirection.create(direction) # to avoid __hash__ calls from a dictionary lookup, # we can inline the small number of options for performance if direction is TransformDirection.FORWARD: return PJ_FWD if direction is TransformDirection.INVERSE: return PJ_INV if direction is TransformDirection.IDENT: return PJ_IDENT raise KeyError(f"{direction} is not a valid TransformDirection") cdef class _TransformerGroup: def __cinit__(self): self.context = NULL self._context_manager = None self._transformers = [] self._unavailable_operations = [] self._best_available = True def __init__( self, _CRS crs_from not None, _CRS crs_to not None, bint always_xy, area_of_interest, bint allow_ballpark, str authority, double accuracy, bint allow_superseded, ): """ From PROJ docs: The operations are sorted with the most relevant ones first: by descending area (intersection of the transformation area with the area of interest, or intersection of the transformation with the area of use of the CRS), and by increasing accuracy. Operations with unknown accuracy are sorted last, whatever their area. """ self.context = pyproj_context_create() self._context_manager = get_context_manager() cdef: PJ_OPERATION_FACTORY_CONTEXT* operation_factory_context = NULL PJ_OBJ_LIST * pj_operations = NULL PJ* pj_transform = NULL const char* c_authority = NULL int num_operations = 0 int is_instantiable = 0 double west_lon_degree double south_lat_degree double east_lon_degree double north_lat_degree if authority is not None: c_authority = authority try: operation_factory_context = proj_create_operation_factory_context( self.context, c_authority, ) if area_of_interest is not None: if not isinstance(area_of_interest, AreaOfInterest): raise ProjError( "Area of interest must be of the type " "pyproj.transformer.AreaOfInterest." ) west_lon_degree = area_of_interest.west_lon_degree south_lat_degree = area_of_interest.south_lat_degree east_lon_degree = area_of_interest.east_lon_degree north_lat_degree = area_of_interest.north_lat_degree proj_operation_factory_context_set_area_of_interest( self.context, operation_factory_context, west_lon_degree, south_lat_degree, east_lon_degree, north_lat_degree, ) if accuracy > 0: proj_operation_factory_context_set_desired_accuracy( self.context, operation_factory_context, accuracy, ) proj_operation_factory_context_set_allow_ballpark_transformations( self.context, operation_factory_context, allow_ballpark, ) proj_operation_factory_context_set_discard_superseded( self.context, operation_factory_context, not allow_superseded, ) proj_operation_factory_context_set_grid_availability_use( self.context, operation_factory_context, PROJ_GRID_AVAILABILITY_IGNORED, ) proj_operation_factory_context_set_spatial_criterion( self.context, operation_factory_context, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION ) pj_operations = proj_create_operations( self.context, crs_from.projobj, crs_to.projobj, operation_factory_context, ) num_operations = proj_list_get_count(pj_operations) for iii in range(num_operations): pj_transform = proj_list_get( self.context, pj_operations, iii, ) is_instantiable = proj_coordoperation_is_instantiable( self.context, pj_transform, ) if is_instantiable: self._transformers.append( _Transformer._from_pj( self.context, pj_transform, always_xy, ) ) else: coordinate_operation = CoordinateOperation.create( self.context, pj_transform, ) self._unavailable_operations.append(coordinate_operation) if iii == 0: self._best_available = False warnings.warn( "Best transformation is not available due to missing " f"{coordinate_operation.grids[0]!r}" ) finally: if operation_factory_context != NULL: proj_operation_factory_context_destroy(operation_factory_context) if pj_operations != NULL: proj_list_destroy(pj_operations) _clear_proj_error() cdef PJ* proj_create_crs_to_crs( PJ_CONTEXT *ctx, const char *source_crs_str, const char *target_crs_str, PJ_AREA *area, str authority, str accuracy, allow_ballpark, bint force_over, only_best, ) except NULL: """ This is the same as proj_create_crs_to_crs in proj.h with the options added. It is a hack for stabilily reasons. Reference: https://github.com/pyproj4/pyproj/pull/800 """ cdef PJ *source_crs = proj_create(ctx, source_crs_str) if source_crs == NULL: _LOGGER.debug( "PROJ_DEBUG: proj_create_crs_to_crs: Cannot instantiate source_crs" ) return NULL cdef PJ *target_crs = proj_create(ctx, target_crs_str) if target_crs == NULL: proj_destroy(source_crs) _LOGGER.debug( "PROJ_DEBUG: proj_create_crs_to_crs: Cannot instantiate target_crs" ) return NULL cdef: const char* options[6] bytes b_authority bytes b_accuracy int options_index = 0 int options_init_iii = 0 for options_init_iii in range(6): options[options_init_iii] = NULL if authority is not None: b_authority = cstrencode(f"AUTHORITY={authority}") options[options_index] = b_authority options_index += 1 if accuracy is not None: b_accuracy = cstrencode(f"ACCURACY={accuracy}") options[options_index] = b_accuracy options_index += 1 if allow_ballpark is not None: if not allow_ballpark: options[options_index] = b"ALLOW_BALLPARK=NO" options_index += 1 if force_over: options[options_index] = b"FORCE_OVER=YES" options_index += 1 if only_best is not None: if only_best: options[options_index] = b"ONLY_BEST=YES" else: options[options_index] = b"ONLY_BEST=NO" cdef PJ* transform = NULL with nogil: transform = proj_create_crs_to_crs_from_pj( ctx, source_crs, target_crs, area, options, ) proj_destroy(source_crs) proj_destroy(target_crs) if transform == NULL: raise ProjError("Error creating Transformer from CRS.") return transform cdef class _Transformer(Base): def __cinit__(self): self._area_of_use = None self.type_name = "Unknown Transformer" self._operations = None self._source_crs = None self._target_crs = None def _initialize_from_projobj(self): self.proj_info = proj_pj_info(self.projobj) if self.proj_info.id == NULL: raise ProjError("Input is not a transformation.") cdef PJ_TYPE transformer_type = proj_get_type(self.projobj) self.type_name = _TRANSFORMER_TYPE_MAP[transformer_type] self._set_base_info() _clear_proj_error() @property def id(self): return self.proj_info.id @property def description(self): return self.proj_info.description @property def definition(self): return self.proj_info.definition @property def has_inverse(self): return self.proj_info.has_inverse == 1 @property def accuracy(self): return self.proj_info.accuracy @property def area_of_use(self): """ Returns ------- AreaOfUse: The area of use object with associated attributes. """ if self._area_of_use is not None: return self._area_of_use self._area_of_use = create_area_of_use(self.context, self.projobj) return self._area_of_use @property def source_crs(self): """ .. versionadded:: 3.3.0 Returns ------- _CRS | None: The source CRS of a CoordinateOperation. """ if self._source_crs is not None: return None if self._source_crs is False else self._source_crs cdef PJ * projobj = proj_get_source_crs(self.context, self.projobj) _clear_proj_error() if projobj == NULL: self._source_crs = False return None try: self._source_crs = _CRS(_to_wkt( self.context, projobj, version=WktVersion.WKT2_2019, pretty=False, )) finally: proj_destroy(projobj) return self._source_crs @property def target_crs(self): """ .. versionadded:: 3.3.0 Returns ------- _CRS | None: The target CRS of a CoordinateOperation. """ if self._target_crs is not None: return None if self._target_crs is False else self._target_crs cdef PJ * projobj = proj_get_target_crs(self.context, self.projobj) _clear_proj_error() if projobj == NULL: self._target_crs = False return None try: self._target_crs = _CRS(_to_wkt( self.context, projobj, version=WktVersion.WKT2_2019, pretty=False, )) finally: proj_destroy(projobj) return self._target_crs @property def operations(self): """ .. versionadded:: 2.4.0 Returns ------- tuple[CoordinateOperation]: The operations in a concatenated operation. """ if self._operations is not None: return self._operations self._operations = _get_concatenated_operations(self.context, self.projobj) return self._operations def get_last_used_operation(self): cdef PJ* last_used_operation = proj_trans_get_last_used_operation(self.projobj) if last_used_operation == NULL: raise ProjError( "Last used operation not found. " "This is likely due to not initiating a transform." ) cdef PJ_CONTEXT* context = NULL try: context = pyproj_context_create() except: proj_destroy(last_used_operation) raise proj_assign_context(last_used_operation, context) return _Transformer._from_pj( context, last_used_operation, False, ) @property def is_network_enabled(self): """ .. versionadded:: 3.0.0 Returns ------- bool: If the network is enabled. """ return proj_context_is_network_enabled(self.context) == 1 def to_proj4(self, version=ProjVersion.PROJ_5, bint pretty=False): """ Convert the projection to a PROJ string. .. versionadded:: 3.1.0 Parameters ---------- version: pyproj.enums.ProjVersion, default=pyproj.enums.ProjVersion.PROJ_5 The version of the PROJ string output. pretty: bool, default=False If True, it will set the output to be a multiline string. Returns ------- str: The PROJ string. """ return _to_proj4(self.context, self.projobj, version=version, pretty=pretty) @staticmethod def from_crs( const char* crs_from, const char* crs_to, bint always_xy=False, area_of_interest=None, str authority=None, str accuracy=None, allow_ballpark=None, bint force_over=False, only_best=None, ): """ Create a transformer from CRS objects """ cdef: PJ_AREA *pj_area_of_interest = NULL double west_lon_degree double south_lat_degree double east_lon_degree double north_lat_degree _Transformer transformer = _Transformer() try: if area_of_interest is not None: if not isinstance(area_of_interest, AreaOfInterest): raise ProjError( "Area of interest must be of the type " "pyproj.transformer.AreaOfInterest." ) pj_area_of_interest = proj_area_create() west_lon_degree = area_of_interest.west_lon_degree south_lat_degree = area_of_interest.south_lat_degree east_lon_degree = area_of_interest.east_lon_degree north_lat_degree = area_of_interest.north_lat_degree proj_area_set_bbox( pj_area_of_interest, west_lon_degree, south_lat_degree, east_lon_degree, north_lat_degree, ) transformer.context = pyproj_context_create() transformer._context_manager = get_context_manager() transformer.projobj = proj_create_crs_to_crs( transformer.context, crs_from, crs_to, pj_area_of_interest, authority=authority, accuracy=accuracy, allow_ballpark=allow_ballpark, force_over=force_over, only_best=only_best, ) finally: if pj_area_of_interest != NULL: proj_area_destroy(pj_area_of_interest) transformer._init_from_crs(always_xy) return transformer @staticmethod cdef _Transformer _from_pj( PJ_CONTEXT* context, PJ *transform_pj, bint always_xy, ): """ Create a Transformer from a PJ* object """ cdef _Transformer transformer = _Transformer() transformer.context = context transformer._context_manager = get_context_manager() transformer.projobj = transform_pj if transformer.projobj == NULL: raise ProjError("Error creating Transformer.") transformer._init_from_crs(always_xy) return transformer @staticmethod def from_pipeline(const char *proj_pipeline): """ Create Transformer from a PROJ pipeline string. """ cdef _Transformer transformer = _Transformer() transformer.context = pyproj_context_create() transformer._context_manager = get_context_manager() auth_match = _AUTH_CODE_RE.match(proj_pipeline.strip()) if auth_match: # attempt to create coordinate operation from AUTH:CODE match_data = auth_match.groupdict() transformer.projobj = proj_create_from_database( transformer.context, cstrencode(match_data["authority"]), cstrencode(match_data["code"]), PJ_CATEGORY_COORDINATE_OPERATION, False, NULL, ) if transformer.projobj == NULL: # initialize projection transformer.projobj = proj_create( transformer.context, proj_pipeline, ) if transformer.projobj is NULL: raise ProjError(f"Invalid projection {proj_pipeline}.") transformer._initialize_from_projobj() return transformer def _set_always_xy(self): """ Setup the transformer so it has the axis order always in xy order. """ cdef PJ* always_xy_pj = proj_normalize_for_visualization( self.context, self.projobj, ) proj_destroy(self.projobj) self.projobj = always_xy_pj def _init_from_crs(self, bint always_xy): """ Finish initializing transformer properties from CRS objects """ if always_xy: self._set_always_xy() self._initialize_from_projobj() @cython.boundscheck(False) @cython.wraparound(False) def _transform( self, object inx, object iny, object inz, object intime, object direction, bint radians, bint errcheck, ): if self.id == "noop": return cdef: PJ_DIRECTION pj_direction = get_pj_direction(direction) PyBuffWriteManager xbuff = PyBuffWriteManager(inx) PyBuffWriteManager ybuff = PyBuffWriteManager(iny) PyBuffWriteManager zbuff PyBuffWriteManager tbuff Py_ssize_t buflenz Py_ssize_t buflent double* zz double* tt if inz is not None: zbuff = PyBuffWriteManager(inz) buflenz = zbuff.len zz = zbuff.data else: buflenz = xbuff.len zz = NULL if intime is not None: tbuff = PyBuffWriteManager(intime) buflent = tbuff.len tt = tbuff.data else: buflent = xbuff.len tt = NULL if not (xbuff.len == ybuff.len == buflenz == buflent): raise ProjError('x, y, z, and time must be same size if included.') cdef Py_ssize_t iii cdef int errno = 0 with nogil: # degrees to radians if not radians and proj_angular_input(self.projobj, pj_direction): for iii in range(xbuff.len): xbuff.data[iii] = xbuff.data[iii]*_DG2RAD ybuff.data[iii] = ybuff.data[iii]*_DG2RAD # radians to degrees elif radians and proj_degree_input(self.projobj, pj_direction): for iii in range(xbuff.len): xbuff.data[iii] = xbuff.data[iii]*_RAD2DG ybuff.data[iii] = ybuff.data[iii]*_RAD2DG proj_errno_reset(self.projobj) proj_trans_generic( self.projobj, pj_direction, xbuff.data, _DOUBLESIZE, xbuff.len, ybuff.data, _DOUBLESIZE, ybuff.len, zz, _DOUBLESIZE, xbuff.len, tt, _DOUBLESIZE, xbuff.len, ) errno = proj_errno(self.projobj) if errcheck and errno: with gil: raise ProjError( f"transform error: {proj_context_errno_string(self.context, errno)}" ) elif errcheck: with gil: if _get_proj_error() is not None: raise ProjError("transform error") # radians to degrees if not radians and proj_angular_output(self.projobj, pj_direction): for iii in range(xbuff.len): xbuff.data[iii] = xbuff.data[iii]*_RAD2DG ybuff.data[iii] = ybuff.data[iii]*_RAD2DG # degrees to radians elif radians and proj_degree_output(self.projobj, pj_direction): for iii in range(xbuff.len): xbuff.data[iii] = xbuff.data[iii]*_DG2RAD ybuff.data[iii] = ybuff.data[iii]*_DG2RAD _clear_proj_error() @cython.boundscheck(False) @cython.wraparound(False) def _transform_point( self, object inx, object iny, object inz, object intime, object direction, bint radians, bint errcheck, ): """ Optimized to transform a single point between two coordinate systems. """ cdef: double coord_x = inx double coord_y = iny double coord_z = 0 double coord_t = HUGE_VAL tuple expected_numeric_types = (int, float) # We do the type-checking internally here due to automatically # casting length-1 arrays to float that we don't want to return scalar for. # Ex: float(np.array([0])) works and we don't want to accept numpy arrays if not isinstance(inx, expected_numeric_types): raise TypeError("Scalar input expected for x") if not isinstance(iny, expected_numeric_types): raise TypeError("Scalar input expected for y") if inz is not None: if not isinstance(inz, expected_numeric_types): raise TypeError("Scalar input expected for z") coord_z = inz if intime is not None: if not isinstance(intime, expected_numeric_types): raise TypeError("Scalar input expected for t") coord_t = intime cdef tuple return_data if self.id == "noop": return_data = (inx, iny) if inz is not None: return_data += (inz,) if intime is not None: return_data += (intime,) return return_data cdef: PJ_DIRECTION pj_direction = get_pj_direction(direction) PJ_COORD projxyout PJ_COORD projxyin = proj_coord(coord_x, coord_y, coord_z, coord_t) with nogil: # degrees to radians if not radians and proj_angular_input(self.projobj, pj_direction): projxyin.uv.u *= _DG2RAD projxyin.uv.v *= _DG2RAD # radians to degrees elif radians and proj_degree_input(self.projobj, pj_direction): projxyin.uv.u *= _RAD2DG projxyin.uv.v *= _RAD2DG proj_errno_reset(self.projobj) projxyout = proj_trans(self.projobj, pj_direction, projxyin) errno = proj_errno(self.projobj) if errcheck and errno: with gil: raise ProjError( f"transform error: {proj_context_errno_string(self.context, errno)}" ) elif errcheck: with gil: if _clear_proj_error() is not None: raise ProjError("transform error") # radians to degrees if not radians and proj_angular_output(self.projobj, pj_direction): projxyout.xy.x *= _RAD2DG projxyout.xy.y *= _RAD2DG # degrees to radians elif radians and proj_degree_output(self.projobj, pj_direction): projxyout.xy.x *= _DG2RAD projxyout.xy.y *= _DG2RAD _clear_proj_error() return_data = (projxyout.xyzt.x, projxyout.xyzt.y) if inz is not None: return_data += (projxyout.xyzt.z,) if intime is not None: return_data += (projxyout.xyzt.t,) return return_data @cython.boundscheck(False) @cython.wraparound(False) def _transform_sequence( self, Py_ssize_t stride, object inseq, bint switch, object direction, bint time_3rd, bint radians, bint errcheck, ): # private function to itransform function if self.id == "noop": return cdef: PJ_DIRECTION pj_direction = get_pj_direction(direction) double *x double *y double *z double *tt if stride < 2: raise ProjError("coordinates must contain at least 2 values") cdef: PyBuffWriteManager coordbuff = PyBuffWriteManager(inseq) Py_ssize_t npts Py_ssize_t iii Py_ssize_t jjj int errno = 0 npts = coordbuff.len // stride with nogil: # degrees to radians if not radians and proj_angular_input(self.projobj, pj_direction): for iii in range(npts): jjj = stride * iii coordbuff.data[jjj] *= _DG2RAD coordbuff.data[jjj + 1] *= _DG2RAD # radians to degrees elif radians and proj_degree_input(self.projobj, pj_direction): for iii in range(npts): jjj = stride * iii coordbuff.data[jjj] *= _RAD2DG coordbuff.data[jjj + 1] *= _RAD2DG if not switch: x = coordbuff.data y = coordbuff.data + 1 else: x = coordbuff.data + 1 y = coordbuff.data # z coordinate if stride == 4 or (stride == 3 and not time_3rd): z = coordbuff.data + 2 else: z = NULL # time if stride == 3 and time_3rd: tt = coordbuff.data + 2 elif stride == 4: tt = coordbuff.data + 3 else: tt = NULL proj_errno_reset(self.projobj) proj_trans_generic( self.projobj, pj_direction, x, stride*_DOUBLESIZE, npts, y, stride*_DOUBLESIZE, npts, z, stride*_DOUBLESIZE, npts, tt, stride*_DOUBLESIZE, npts, ) errno = proj_errno(self.projobj) if errcheck and errno: with gil: raise ProjError( f"itransform error: {proj_context_errno_string(self.context, errno)}" ) elif errcheck: with gil: if _get_proj_error() is not None: raise ProjError("itransform error") # radians to degrees if not radians and proj_angular_output(self.projobj, pj_direction): for iii in range(npts): jjj = stride * iii coordbuff.data[jjj] *= _RAD2DG coordbuff.data[jjj + 1] *= _RAD2DG # degrees to radians elif radians and proj_degree_output(self.projobj, pj_direction): for iii in range(npts): jjj = stride * iii coordbuff.data[jjj] *= _DG2RAD coordbuff.data[jjj + 1] *= _DG2RAD _clear_proj_error() @cython.boundscheck(False) @cython.wraparound(False) def _transform_bounds( self, double left, double bottom, double right, double top, int densify_pts, bint radians, bint errcheck, object direction, ): cdef PJ_DIRECTION pj_direction = get_pj_direction(direction) if self.id == "noop" or pj_direction == PJ_IDENT: return (left, bottom, right, top) cdef: int errno = 0 bint success = True double out_left = left double out_bottom = bottom double out_right = right double out_top = top with nogil: # degrees to radians if not radians and proj_angular_input(self.projobj, pj_direction): left *= _DG2RAD bottom *= _DG2RAD right *= _DG2RAD top *= _DG2RAD # radians to degrees elif radians and proj_degree_input(self.projobj, pj_direction): left *= _RAD2DG bottom *= _RAD2DG right *= _RAD2DG top *= _RAD2DG proj_errno_reset(self.projobj) success = proj_trans_bounds( self.context, self.projobj, pj_direction, left, bottom, right, top, &out_left, &out_bottom, &out_right, &out_top, densify_pts, ) if not success or errcheck: errno = proj_errno(self.projobj) if errno: with gil: raise ProjError( "transform bounds error: " f"{proj_context_errno_string(self.context, errno)}" ) else: with gil: if _get_proj_error() is not None: raise ProjError("transform bounds error") # radians to degrees if not radians and proj_angular_output(self.projobj, pj_direction): out_left *= _RAD2DG out_bottom *= _RAD2DG out_right *= _RAD2DG out_top *= _RAD2DG # degrees to radians elif radians and proj_degree_output(self.projobj, pj_direction): out_left *= _DG2RAD out_bottom *= _DG2RAD out_right *= _DG2RAD out_top *= _DG2RAD _clear_proj_error() return out_left, out_bottom, out_right, out_top @cython.boundscheck(False) @cython.wraparound(False) def _get_factors(self, longitude, latitude, bint radians, bint errcheck): """ Calculates the projection factors PJ_FACTORS Designed to work with Proj class. Equivalent to `proj -S` command line. """ cdef PyBuffWriteManager lonbuff = PyBuffWriteManager(longitude) cdef PyBuffWriteManager latbuff = PyBuffWriteManager(latitude) if not lonbuff.len or not (lonbuff.len == latbuff.len): raise ProjError('longitude and latitude must be same size') # prepare the factors output meridional_scale = copy.copy(longitude) parallel_scale = copy.copy(longitude) areal_scale = copy.copy(longitude) angular_distortion = copy.copy(longitude) meridian_parallel_angle = copy.copy(longitude) meridian_convergence = copy.copy(longitude) tissot_semimajor = copy.copy(longitude) tissot_semiminor = copy.copy(longitude) dx_dlam = copy.copy(longitude) dx_dphi = copy.copy(longitude) dy_dlam = copy.copy(longitude) dy_dphi = copy.copy(longitude) cdef: PyBuffWriteManager meridional_scale_buff = PyBuffWriteManager( meridional_scale ) PyBuffWriteManager parallel_scale_buff = PyBuffWriteManager( parallel_scale ) PyBuffWriteManager areal_scale_buff = PyBuffWriteManager(areal_scale) PyBuffWriteManager angular_distortion_buff = PyBuffWriteManager( angular_distortion ) PyBuffWriteManager meridian_parallel_angle_buff = PyBuffWriteManager( meridian_parallel_angle ) PyBuffWriteManager meridian_convergence_buff = PyBuffWriteManager( meridian_convergence ) PyBuffWriteManager tissot_semimajor_buff = PyBuffWriteManager( tissot_semimajor ) PyBuffWriteManager tissot_semiminor_buff = PyBuffWriteManager( tissot_semiminor ) PyBuffWriteManager dx_dlam_buff = PyBuffWriteManager(dx_dlam) PyBuffWriteManager dx_dphi_buff = PyBuffWriteManager(dx_dphi) PyBuffWriteManager dy_dlam_buff = PyBuffWriteManager(dy_dlam) PyBuffWriteManager dy_dphi_buff = PyBuffWriteManager(dy_dphi) # calculate the factors PJ_COORD pj_coord = proj_coord(0, 0, 0, HUGE_VAL) PJ_FACTORS pj_factors int errno = 0 bint invalid_coord = 0 Py_ssize_t iii with nogil: for iii in range(lonbuff.len): pj_coord.uv.u = lonbuff.data[iii] pj_coord.uv.v = latbuff.data[iii] if not radians: pj_coord.uv.u *= _DG2RAD pj_coord.uv.v *= _DG2RAD # set both to HUGE_VAL if inf or nan proj_errno_reset(self.projobj) if pj_coord.uv.v == HUGE_VAL \ or pj_coord.uv.v != pj_coord.uv.v \ or pj_coord.uv.u == HUGE_VAL \ or pj_coord.uv.u != pj_coord.uv.u: invalid_coord = True else: invalid_coord = False pj_factors = proj_factors(self.projobj, pj_coord) errno = proj_errno(self.projobj) if errcheck and errno: with gil: raise ProjError( f"proj error: {proj_context_errno_string(self.context, errno)}" ) if errno or invalid_coord: meridional_scale_buff.data[iii] = HUGE_VAL parallel_scale_buff.data[iii] = HUGE_VAL areal_scale_buff.data[iii] = HUGE_VAL angular_distortion_buff.data[iii] = HUGE_VAL meridian_parallel_angle_buff.data[iii] = HUGE_VAL meridian_convergence_buff.data[iii] = HUGE_VAL tissot_semimajor_buff.data[iii] = HUGE_VAL tissot_semiminor_buff.data[iii] = HUGE_VAL dx_dlam_buff.data[iii] = HUGE_VAL dx_dphi_buff.data[iii] = HUGE_VAL dy_dlam_buff.data[iii] = HUGE_VAL dy_dphi_buff.data[iii] = HUGE_VAL else: meridional_scale_buff.data[iii] = pj_factors.meridional_scale parallel_scale_buff.data[iii] = pj_factors.parallel_scale areal_scale_buff.data[iii] = pj_factors.areal_scale angular_distortion_buff.data[iii] = ( pj_factors.angular_distortion * _RAD2DG ) meridian_parallel_angle_buff.data[iii] = ( pj_factors.meridian_parallel_angle * _RAD2DG ) meridian_convergence_buff.data[iii] = ( pj_factors.meridian_convergence * _RAD2DG ) tissot_semimajor_buff.data[iii] = pj_factors.tissot_semimajor tissot_semiminor_buff.data[iii] = pj_factors.tissot_semiminor dx_dlam_buff.data[iii] = pj_factors.dx_dlam dx_dphi_buff.data[iii] = pj_factors.dx_dphi dy_dlam_buff.data[iii] = pj_factors.dy_dlam dy_dphi_buff.data[iii] = pj_factors.dy_dphi _clear_proj_error() return Factors( meridional_scale=meridional_scale, parallel_scale=parallel_scale, areal_scale=areal_scale, angular_distortion=angular_distortion, meridian_parallel_angle=meridian_parallel_angle, meridian_convergence=meridian_convergence, tissot_semimajor=tissot_semimajor, tissot_semiminor=tissot_semiminor, dx_dlam=dx_dlam, dx_dphi=dx_dphi, dy_dlam=dy_dlam, dy_dphi=dy_dphi, )