diff --git a/scipy/cluster/_hierarchy.pyx b/scipy/cluster/_hierarchy.pyx
index 373f94a494b2..828344283af0 100644
--- a/scipy/cluster/_hierarchy.pyx
+++ b/scipy/cluster/_hierarchy.pyx
@@ -11,6 +11,7 @@ cdef extern from "numpy/npy_math.h":
 
 ctypedef unsigned char uchar
 
+np.import_array()
 
 # _hierarchy_distance_update.pxi includes the definition of linkage_distance_update
 # and the distance update functions for the supported linkage methods.
diff --git a/scipy/cluster/_optimal_leaf_ordering.pyx b/scipy/cluster/_optimal_leaf_ordering.pyx
index f9b320593fdf..4feaadbc4e73 100644
--- a/scipy/cluster/_optimal_leaf_ordering.pyx
+++ b/scipy/cluster/_optimal_leaf_ordering.pyx
@@ -32,6 +32,7 @@ from libc.stdlib cimport malloc, free
 
 from scipy.spatial.distance import squareform, is_valid_y, is_valid_dm
 
+np.import_array()
 
 @cython.profile(False)
 @cython.boundscheck(False)
diff --git a/scipy/fftpack/convolve.pyx b/scipy/fftpack/convolve.pyx
index 2dc85c356065..76a58f2dcf0e 100644
--- a/scipy/fftpack/convolve.pyx
+++ b/scipy/fftpack/convolve.pyx
@@ -3,6 +3,8 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
+np.import_array()
+
 __all__ = ['destroy_convolve_cache', 'convolve', 'convolve_z',
            'init_convolution_kernel']
 
diff --git a/scipy/interpolate/_bspl.pyx b/scipy/interpolate/_bspl.pyx
index 5097b58b265b..f1c0b126110c 100644
--- a/scipy/interpolate/_bspl.pyx
+++ b/scipy/interpolate/_bspl.pyx
@@ -8,6 +8,8 @@ cimport numpy as cnp
 
 cimport cython
 
+cnp.import_array()
+
 cdef extern from "src/__fitpack.h":
     void _deBoor_D(const double *t, double x, int k, int ell, int m, double *result) nogil
 
diff --git a/scipy/io/matlab/mio_utils.pyx b/scipy/io/matlab/mio_utils.pyx
index b1ff2788d5bc..56f5c176af25 100644
--- a/scipy/io/matlab/mio_utils.pyx
+++ b/scipy/io/matlab/mio_utils.pyx
@@ -5,6 +5,8 @@
 import numpy as np
 cimport numpy as cnp
 
+cnp.import_array()
+
 
 cpdef object squeeze_element(cnp.ndarray arr):
     ''' Return squeezed element
diff --git a/scipy/ndimage/src/_cytest.pyx b/scipy/ndimage/src/_cytest.pyx
index 44f5ab8178b2..3948d97a78c7 100644
--- a/scipy/ndimage/src/_cytest.pyx
+++ b/scipy/ndimage/src/_cytest.pyx
@@ -6,6 +6,7 @@ from cpython.pycapsule cimport (
 cimport numpy as np
 from numpy cimport npy_intp as intp
 
+np.import_array()
 
 cdef void _destructor(obj):
     cdef void *callback_data = PyCapsule_GetContext(obj)
diff --git a/scipy/optimize/_bglu_dense.pyx b/scipy/optimize/_bglu_dense.pyx
index 8ae4889a2a75..376df50ad125 100644
--- a/scipy/optimize/_bglu_dense.pyx
+++ b/scipy/optimize/_bglu_dense.pyx
@@ -11,6 +11,8 @@ try:
 except ImportError:
     from time import clock as timer
 
+np.import_array()
+
 __all__ = ['LU', 'BGLU']
 
 @cython.boundscheck(False)
diff --git a/scipy/optimize/_group_columns.pyx b/scipy/optimize/_group_columns.pyx
index 7b6e933d4639..510315c9201f 100644
--- a/scipy/optimize/_group_columns.pyx
+++ b/scipy/optimize/_group_columns.pyx
@@ -10,6 +10,8 @@ import numpy as np
 cimport numpy as np
 from cpython cimport bool
 
+np.import_array()
+
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
diff --git a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
index 579f7e057ac1..825c77b2e233 100644
--- a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
+++ b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx
@@ -64,6 +64,8 @@ from .HighsOptions cimport (
 )
 from .HighsModelUtils cimport utilPrimalDualStatusToString
 
+np.import_array()
+
 # options to reference for default values and bounds;
 # make a map to quickly lookup
 cdef HighsOptions _ref_opts
diff --git a/scipy/optimize/_trlib/_trlib.pyx b/scipy/optimize/_trlib/_trlib.pyx
index 1f1a513ba5b5..9aa9cd7a2f03 100644
--- a/scipy/optimize/_trlib/_trlib.pyx
+++ b/scipy/optimize/_trlib/_trlib.pyx
@@ -6,6 +6,7 @@ cimport numpy as np
 
 from scipy._lib.messagestream cimport MessageStream
 
+np.import_array()
 
 class TRLIBQuadraticSubproblem(BaseQuadraticSubproblem):
 
diff --git a/scipy/signal/_max_len_seq_inner.pyx b/scipy/signal/_max_len_seq_inner.pyx
index 7f448a1d8024..8ce5bc4e72a2 100644
--- a/scipy/signal/_max_len_seq_inner.pyx
+++ b/scipy/signal/_max_len_seq_inner.pyx
@@ -5,6 +5,7 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
+np.import_array()
 
 # Fast inner loop of max_len_seq.
 @cython.cdivision(True)  # faster modulo
diff --git a/scipy/signal/_peak_finding_utils.pyx b/scipy/signal/_peak_finding_utils.pyx
index 309d72a0ae0e..4a523f135ca2 100644
--- a/scipy/signal/_peak_finding_utils.pyx
+++ b/scipy/signal/_peak_finding_utils.pyx
@@ -11,6 +11,7 @@ import numpy as np
 cimport numpy as np
 from libc.math cimport ceil
 
+np.import_array()
 
 __all__ = ['_local_maxima_1d', '_select_by_peak_distance', '_peak_prominences',
            '_peak_widths']
diff --git a/scipy/signal/_sosfilt.pyx b/scipy/signal/_sosfilt.pyx
index cf626bab29f4..b62070b7dfd2 100644
--- a/scipy/signal/_sosfilt.pyx
+++ b/scipy/signal/_sosfilt.pyx
@@ -1,6 +1,8 @@
 cimport numpy as np
 cimport cython
 
+np.import_array()
+
 ctypedef fused DTYPE_floating_t:
     float
     float complex
diff --git a/scipy/signal/_spectral.pyx b/scipy/signal/_spectral.pyx
index 374a98575a5c..e13c3a96d4d0 100644
--- a/scipy/signal/_spectral.pyx
+++ b/scipy/signal/_spectral.pyx
@@ -7,6 +7,8 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
+np.import_array()
+
 __all__ = ['_lombscargle']
 
 
diff --git a/scipy/signal/_upfirdn_apply.pyx b/scipy/signal/_upfirdn_apply.pyx
index 42fe67921584..0e23594023d0 100644
--- a/scipy/signal/_upfirdn_apply.pyx
+++ b/scipy/signal/_upfirdn_apply.pyx
@@ -40,6 +40,7 @@ from cython import bint  # boolean integer type
 from libc.stdlib cimport malloc, free
 from libc.string cimport memset
 
+np.import_array()
 
 ctypedef double complex double_complex
 ctypedef float complex float_complex
diff --git a/scipy/sparse/csgraph/_flow.pyx b/scipy/sparse/csgraph/_flow.pyx
index 6153a946c6a8..0a2eb4fad77c 100644
--- a/scipy/sparse/csgraph/_flow.pyx
+++ b/scipy/sparse/csgraph/_flow.pyx
@@ -7,6 +7,8 @@ cimport numpy as np
 
 include 'parameters.pxi'
 
+np.import_array()
+
 
 class MaximumFlowResult:
     """Represents the result of a maximum flow calculation.
diff --git a/scipy/sparse/csgraph/_matching.pyx b/scipy/sparse/csgraph/_matching.pyx
index 0f6bcfaa641e..1017b523d830 100644
--- a/scipy/sparse/csgraph/_matching.pyx
+++ b/scipy/sparse/csgraph/_matching.pyx
@@ -9,6 +9,8 @@ from numpy.math cimport INFINITY
 from scipy.sparse import (csr_matrix,
                           isspmatrix_coo, isspmatrix_csc, isspmatrix_csr)
 
+np.import_array()
+
 include "parameters.pxi"
 
 
diff --git a/scipy/sparse/csgraph/_min_spanning_tree.pyx b/scipy/sparse/csgraph/_min_spanning_tree.pyx
index ecdefdfd7b58..2d3f18921332 100644
--- a/scipy/sparse/csgraph/_min_spanning_tree.pyx
+++ b/scipy/sparse/csgraph/_min_spanning_tree.pyx
@@ -8,6 +8,8 @@ cimport cython
 from scipy.sparse import csr_matrix, isspmatrix_csc, isspmatrix
 from scipy.sparse.csgraph._validation import validate_graph
 
+np.import_array()
+
 include 'parameters.pxi'
 
 def minimum_spanning_tree(csgraph, overwrite=False):
diff --git a/scipy/sparse/csgraph/_reordering.pyx b/scipy/sparse/csgraph/_reordering.pyx
index 2223b37827f7..a95d4f1227ea 100644
--- a/scipy/sparse/csgraph/_reordering.pyx
+++ b/scipy/sparse/csgraph/_reordering.pyx
@@ -10,6 +10,7 @@ from scipy.sparse import (csc_matrix, csr_matrix, isspmatrix, isspmatrix_coo,
                           SparseEfficiencyWarning)
 from . import maximum_bipartite_matching
 
+np.import_array()
 
 include 'parameters.pxi'
 
diff --git a/scipy/sparse/csgraph/_shortest_path.pyx b/scipy/sparse/csgraph/_shortest_path.pyx
index 8f1056239e07..97208d055cad 100644
--- a/scipy/sparse/csgraph/_shortest_path.pyx
+++ b/scipy/sparse/csgraph/_shortest_path.pyx
@@ -22,6 +22,8 @@ cimport cython
 from libc.stdlib cimport malloc, free
 from numpy.math cimport INFINITY
 
+np.import_array()
+
 include 'parameters.pxi'
 
 
diff --git a/scipy/sparse/csgraph/_tools.pyx b/scipy/sparse/csgraph/_tools.pyx
index 3ae4bbfa3719..dd2208903a3d 100644
--- a/scipy/sparse/csgraph/_tools.pyx
+++ b/scipy/sparse/csgraph/_tools.pyx
@@ -11,6 +11,8 @@ cimport numpy as np
 from scipy.sparse import csr_matrix, isspmatrix,\
     isspmatrix_csr, isspmatrix_csc, isspmatrix_lil
 
+np.import_array()
+
 include 'parameters.pxi'
 
 def csgraph_from_masked(graph):
diff --git a/scipy/sparse/csgraph/_traversal.pyx b/scipy/sparse/csgraph/_traversal.pyx
index 9261e968fc2a..52ba0de6442f 100644
--- a/scipy/sparse/csgraph/_traversal.pyx
+++ b/scipy/sparse/csgraph/_traversal.pyx
@@ -15,6 +15,8 @@ from scipy.sparse.csgraph._tools import reconstruct_path
 cimport cython
 from libc cimport stdlib
 
+np.import_array()
+
 include 'parameters.pxi'
 
 def connected_components(csgraph, directed=True, connection='weak',
diff --git a/scipy/spatial/_hausdorff.pyx b/scipy/spatial/_hausdorff.pyx
index 0f023d143184..0b96fdf46ff4 100644
--- a/scipy/spatial/_hausdorff.pyx
+++ b/scipy/spatial/_hausdorff.pyx
@@ -15,6 +15,8 @@ cimport numpy as np
 cimport cython
 from libc.math cimport sqrt
 
+np.import_array()
+
 __all__ = ['directed_hausdorff']
 
 @cython.boundscheck(False)
diff --git a/scipy/spatial/_voronoi.pyx b/scipy/spatial/_voronoi.pyx
index e6c09562a88b..f86086a4aac1 100644
--- a/scipy/spatial/_voronoi.pyx
+++ b/scipy/spatial/_voronoi.pyx
@@ -14,6 +14,8 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
+np.import_array()
+
 __all__ = ['sort_vertices_of_regions']
 
 # array-filling placeholder that can never occur
diff --git a/scipy/spatial/ckdtree.pyx b/scipy/spatial/ckdtree.pyx
index 18d3e806e3a7..6de28e191748 100644
--- a/scipy/spatial/ckdtree.pyx
+++ b/scipy/spatial/ckdtree.pyx
@@ -23,6 +23,8 @@ import threading
 import operator
 import warnings
 
+np.import_array()
+
 cdef extern from "<limits.h>":
     long LONG_MAX
 
diff --git a/scipy/spatial/qhull.pyx b/scipy/spatial/qhull.pyx
index 90f207a8821f..828fa21b1913 100644
--- a/scipy/spatial/qhull.pyx
+++ b/scipy/spatial/qhull.pyx
@@ -23,6 +23,8 @@ import os
 import sys
 import tempfile
 
+np.import_array()
+
 cdef extern from "numpy/npy_math.h":
     double nan "NPY_NAN"
 
diff --git a/scipy/spatial/transform/rotation.pyx b/scipy/spatial/transform/rotation.pyx
index 2e5a90d1c775..bc1eba640137 100644
--- a/scipy/spatial/transform/rotation.pyx
+++ b/scipy/spatial/transform/rotation.pyx
@@ -10,6 +10,8 @@ from cython.view cimport array
 from libc.math cimport sqrt, sin, cos, atan2, acos
 from numpy.math cimport PI as pi, NAN, isnan # avoid MSVC error
 
+np.import_array()
+
 # utilities for empty array initialization
 cdef inline double[:] _empty1(int n):
     return array(shape=(n,), itemsize=sizeof(double), format=b"d")
diff --git a/scipy/stats/_sobol.pyx b/scipy/stats/_sobol.pyx
index 702fc95aed13..49c2e154af79 100644
--- a/scipy/stats/_sobol.pyx
+++ b/scipy/stats/_sobol.pyx
@@ -6,6 +6,8 @@ cimport numpy as cnp
 import os
 import numpy as np
 
+cnp.import_array()
+
 # Parameters are linked to the direction numbers list.
 # See `initialize_direction_numbers` for more details.
 # Declared using DEF to be known at compilation time for ``poly`` et ``vinit``
diff --git a/scipy/stats/_stats.pyx b/scipy/stats/_stats.pyx
index 528531e7a22d..7e9744e3b5d1 100644
--- a/scipy/stats/_stats.pyx
+++ b/scipy/stats/_stats.pyx
@@ -12,6 +12,7 @@ import numpy as np
 import scipy.stats, scipy.special
 cimport scipy.special.cython_special as cs
 
+np.import_array()
 
 cdef double von_mises_cdf_series(double k, double x, unsigned int p):
     cdef double s, c, sn, cn, R, V
diff --git a/scipy/stats/biasedurn.pyx.templ b/scipy/stats/biasedurn.pyx.templ
index dfc9ac767f8e..59bc8625c8c7 100644
--- a/scipy/stats/biasedurn.pyx.templ
+++ b/scipy/stats/biasedurn.pyx.templ
@@ -7,6 +7,8 @@ cimport numpy as np
 import numpy as np
 from libcpp.memory cimport unique_ptr
 
+np.import_array()
+
 IF not NPY_OLD:
     from numpy.random cimport bitgen_t
     from cpython.pycapsule cimport PyCapsule_GetPointer, PyCapsule_IsValid