pytroll / pyresample Goto Github PK

Geospatial image resampling in Python

Home Page: http://pyresample.readthedocs.org

License: GNU Lesser General Public License v3.0

Python 97.98% C++ 1.51% CSS 0.22% HTML 0.29%

python numpy resampling kd-tree hacktoberfest dask xarray closember

pyresample's Introduction

This is the sandbox area for the pytroll project, an international cooperation 
on a future distributed real-time processing system for Meteorological Satellite
Data.



------------------------------------
December 2010

Lars Ørum Rasmussen, Esben Sigård Nielsen, Kristian Rune Larsen, 
Martin Raspaud, Anna Geidne, Adam Dybbroe.

Danish Meteorological Institute (DMI)
Swedish Meteorological and Hydrological Institute (SMHI)

pyresample's People

Stargazers

Watchers

Forkers

valgur spareeth djhoese solab jsolbrig mindyls npiaget bhawkins meteoswiss-mdr cpaulik avalentino deni90 shiweihappy moondav brendancol mitkin sunshine1204 khunger aprilwaterlily joycezw goodsonr donkpixi howff geoguxm agoodm exwexs art816 hyoklee sebastic nojuman frdcms anirudh-ramesh wroberts4 drnextgis bertilhatt wk1984 mraspaud pnuu loreclem gerritholl observingclouds leosiqueira adybbroe vikashkodati wsf1990 arpm cycle13 ninahakansson mostamndi kuonanhong anuragsinghchaudhary hundahl ibrewster fagan2888 ifenty tomlav ameraner tsukada-cs ouyangradi tanjased benr0 jleinonen cgentemann nourou6 bihtert zxdawn abrammer deweatherman owenlittlejohns chanjeunlam maxmarumbwa danschef olivier-membrive yunjunz ghiggi simonrp84 bochuxt isabella232 benjaminesse plantain bzah sfinkens flamingbear gschen23 swamydev iq-scm alexandremary ethantkoenig jgirouard23 cmo84

pyresample's Issues

Resampling of swath data - result shape difference

Hello, I'm new to this, but we are looking at using pyresample to speed up our mapping routines of h5 L2 satellite data from polar orbiters. We are currently using GPT provided in BEAM which is java based and very VERY slow.

However, when I tried to use your samples as in the docs:

#!/usr/bin/env python3

import pyresample as pr
from pyresample import kd_tree, image, bilinear, geometry
import numpy as np

area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',{'a': '6378144.0', 'b': '6356759.0','lat_0': '50.00', 'lat_ts': '50.00','lon_0': '8.00', 'proj': 'stere'},800, 800,[-1370912.72, -909968.64,1029087.28, 1490031.36])
data = np.fromfunction(lambda y, x: y*x, (50, 10))
lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)

#using pyresample.image
swath_con = image.ImageContainerNearest(data, swath_def, radius_of_influence=5000)
area_con = swath_con.resample(area_def)
swath_con_result = area_con.image_data

#using pyresample.kd_tree
resample_nearest_result = kd_tree.resample_nearest(swath_def, data, area_def, radius_of_influence=50000,epsilon=10)
resample_gauss_result = kd_tree.resample_gauss(swath_def, data, area_def, radius_of_influence=50000,sigmas=100)

#using pyresample.bilinear (it appears that the gridding does not take place?
resample_bilinear_result = bilinear.resample_bilinear(data, swath_def, area_def,radius=30000, neighbours=15,nprocs=1, fill_value=0,reduce_data=True, segments=None,epsilon=0)

print("Starting data:")
print("data.min(): ", data.min())
print("data.max(): ", data.max())
print()
print("image.ImageContainerNearest method:")
print("swath_con_result.shape: ", swath_con_result.shape)
print("swath_con_result.size: ", swath_con_result.size)
print("swath_con_result.min(): ",swath_con_result.min())
print("swath_con_result.max(): ",swath_con_result.max())
print()
print("kd_tree.resample_nearest method:")
print("resample_nearest_result.shape: ", resample_nearest_result.shape)
print("resample_nearest_result.size: ", resample_nearest_result.size)
print("resample_nearest_result.min(): ", resample_nearest_result.min())
print("resample_nearest_result.max(): ", resample_nearest_result.max())
print()
print("kd_tree.resample_gauss method:")
print("resample_gauss_result.shape: ", resample_gauss_result.shape)
print("resample_gauss_result.size: ", resample_gauss_result.size)
print("resample_gauss_result.min(): ", resample_gauss_result.min())
print("resample_gauss_result.max(): ", resample_gauss_result.max())
print()
print("bilinear.resample_bilinear method:")
print("resample_bilinear_result.shape: ", resample_bilinear_result.shape)
print("resample_bilinear_result.size: ", resample_bilinear_result.size)
print("resample_bilinear_result.min(): ", resample_bilinear_result.min())
print("resample_bilinear_result.max(): ", resample_bilinear_result.max())

The results.shape is unexpected for the resample_bilinear_result when I run the resample_test.py file from the code above eg:

./resample_tests.py
Starting data:
data.min(): 0.0
data.max(): 441.0

image.ImageContainerNearest method:
swath_con_result.shape: (800, 800)
swath_con_result.size: 640000
swath_con_result.min(): 0.0
swath_con_result.max(): 297.0

kd_tree.resample_nearest method:
resample_nearest_result.shape: (800, 800)
resample_nearest_result.size: 640000
resample_nearest_result.min(): 0.0
resample_nearest_result.max(): 297.0

kd_tree.resample_gauss method:
resample_gauss_result.shape: (800, 800)
resample_gauss_result.size: 640000
resample_gauss_result.min(): 0.0
resample_gauss_result.max(): 297.0

bilinear.resample_bilinear method:
resample_bilinear_result.shape: (640000, 1)
resample_bilinear_result.size: 640000
resample_bilinear_result.min(): 0.0
resample_bilinear_result.max(): 0.0

And regardless of what I change in the bilinear.resample_bilinear method, I always get 0 min and max vals and I'm only guessing, but is that because the data isn't gridded at this point?

Thanks,
Brock

test_files not included in the source package

test_files not included in the source package on PyPI.
Without test files it is not possible to run the test suite.
Please add missing files or remove the test suite form the source distribution.

sys.path changed during install

From jakob.malm on August 03, 2011 17:03:15

I would like to install pyresample automatically as a dependency in my setuptools setup.py script. setup.py of pyresample unfortunately inserts a path at the first position in sys.path, such that any subsequent imports of 'version' imports the one from pyresample.

Please reset sys.path after use.

Original issue: http://code.google.com/p/pyresample/issues/detail?id=3

Allow proj strings for area definition intialization

Currently AreaDefinition accepts only proj4 dicts. Would be useful to pass proj4 strings, since Proj can parse both.

Versions from 1.2.0 fail to build on Windows.

I'm no expert on Windows compilation but noticed in CI that pyresample fails to build for versions greater than 1.1.6.

See here for a successful build and here for a broken one.

You might want to consider using Appveyor to test your software on Windows if you want to support it.

Wrong formula used for the full width at half maximum

In pyresample/utils.py, the formula in the function fwhm2sigma() is as follow:
fwhm / (2 * np.sqrt(np.log(2))).
However, it should be:
fwhm / (2 * np.sqrt(2 * np.log(2))),
according to Mathworld: Gaussian Function Equation 8.

kd_tree.get_neighbour_info won't accept fill_value keyword.

It is more of a documentation issue than the code. It looks like kd_tree.get_neighbour_info can't accept fill_value keyword despite the docstring saying so.

fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked

I guess the easiest fix is to simply edit the docstring.

ImportError: undefined symbol: PyFPE_jbuf

Filing an issue for the reference:

I experienced an error when importing fornav and ll2cr:

    importError: /home/mikhail/.virtualenvs/def27/local/lib/python2.7/site-packages/pyresample/ewa/_ll2cr.so: undefined symbol: PyFPE_jbuf

Apparently this is caused by two incompatible versions of numpy installed, one in virtualenv and another as a system package python-numpy. Obviously the virtualenv was not totally isolated from system site-packages.

ImageContainerQuick does not mask fill values in high resolution images

ImageContainerQuick.resample() does not mask fill values in high resolution images. The output is a masked array, but its mask is False. It works fine with low resolutions though. The problem seems to be independent of the projection (I didn't check all).

I'm using pyresample-1.5.0 with numpy-1.13.0 .

Example:

import numpy as np
import pyresample.geometry
import pyresample.image
import pyresample.plot
import matplotlib.pyplot as plt
from pycmsaf.geometry import get_proj_corners

for nx, ny in zip([360, 3600], [180, 1800]):
    # Generate some test data with masked elements
    x = np.arange(nx, dtype='f4')
    y = np.arange(ny, dtype='f4')
    x *= 2 * np.pi / x.max()
    y *= 2 * np.pi / y.max()
    mx, my = np.meshgrid(x, y)
    data = np.sin(np.sqrt(mx ** 2 + my ** 2))
    data = np.ma.masked_where(data > 0.98, data)

    # Define regular lon-lat grid (lower left corner: [-180, -90], upper right
    # corner: [180, 90])
    regular_grid = pyresample.geometry.AreaDefinition(
        area_id='my_regular_grid',
        name='Regular 1x1 degree grid',
        proj_id='my_id',
        proj_dict={'proj': 'eqc', 'lon_0': 0},
        x_size=nx,
        y_size=ny,
        area_extent=[-20037508.342789244, -10018754.171394622,
                     20037508.342789244, 10018754.171394622])

    # Define web-mercator grid (lower left corner: [-180, -85], upper right
    # corner: [180, 85])
    proj4_str = '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 ' \
                '+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null'
    web_mercator = pyresample.utils.get_area_def(
        area_id='my_web_mercator',
        area_name='Web Mercator',
        proj_id='my_id',
        proj4_args=proj4_str,
        x_size=ny,
        y_size=ny,
        area_extent=[-20037508.342789244, -19971868.88040857,
                     20037508.342789244, 19971868.880408563]
    )

    # Resample data from regular lon-lat grid to web-mercator grid
    con_reg = pyresample.image.ImageContainerQuick(
        image_data=data, geo_def=regular_grid, fill_value=None)
    con_merc = con_reg.resample(web_mercator)
    data_merc = con_merc.image_data
    print data_merc.mask

    # Plot data
    fig, ax = plt.subplots()
    bmap = pyresample.plot.area_def2basemap(web_mercator, resolution='c')
    bmap.ax = ax
    bmap.drawcoastlines()
    bmap.imshow(data_merc, origin='upper')
    ax.set_title('nx={0}, ny={1}'.format(nx, ny))

plt.show()

Output:

Let the user provide custom weight into pyresample

From [email protected] on November 25, 2013 14:15:49

Hei pyresample team,

I am seating at the pytroll workshop, and would like a new feature:

Users could specify an array of weights (same dimension as lat/lon/data), and this weight would multiply the weight computed from distance (through the fov function, either _custom() or _gauss()).

Two use cases:

the weight could be given as 1./var(data) : the inverse of the uncertainty in each data point (as variance);
the weight could be computed on fov observation time, with more weight towards 12utc than at 0utc and 23:59utc. A resulting daily composite image would then be "sharpened towards 12utc", wrt to a daily composite without any weighting.

I'll start implementing right now :)
Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=12

Possible to remap from grid to swath?

From jakob.malm on May 19, 2010 16:56:31

I had a quick look in the documentation but did not find a way to remap
gridded data to swath. Is this possible? I would like to remap NWP data to
a satellite overpass (with full geolocation) in an easy way.

Original issue: http://code.google.com/p/pyresample/issues/detail?id=1

Test failure on debian sid x86_64

I'm try to generate a new debian package for pyresample.

$ uname -a
Linux mac2 3.19.0-23-generic #24-Ubuntu SMP Tue Jul 7 18:52:55 UTC 2015 x86_64 GNU/Linux

When I run the test suite I get 5 errors and some deprecation warning:

$ env PYTHONPATH=. python2.7 /usr/bin/nosetests -v -w pyresample/test

test_append_1d (pyresample.test.test_geometry.Test) ... ok
test_append_2d (pyresample.test.test_geometry.Test) ... ok
test_area_equal (pyresample.test.test_geometry.Test) ... ok
test_area_extent_ll (pyresample.test.test_geometry.Test) ... ok
test_base_lat_invalid (pyresample.test.test_geometry.Test) ... ok
test_base_lon_wrapping (pyresample.test.test_geometry.Test) ... ok
test_base_type (pyresample.test.test_geometry.Test) ... ok
test_boundary (pyresample.test.test_geometry.Test) ... ok
test_cartesian (pyresample.test.test_geometry.Test) ... ok
test_concat_1d (pyresample.test.test_geometry.Test) ... ok
test_concat_2d (pyresample.test.test_geometry.Test) ... ok
Test the function get_xy_from_lonlat ... ok
test_grid_filter (pyresample.test.test_geometry.Test) ... ERROR
test_grid_filter2D (pyresample.test.test_geometry.Test) ... ERROR
test_grid_filter_valid (pyresample.test.test_geometry.Test) ... ERROR
test_latlong_area (pyresample.test.test_geometry.Test) ... ERROR
test_lonlat_precomp (pyresample.test.test_geometry.Test) ... ok
test_not_area_equal (pyresample.test.test_geometry.Test) ... ok
test_swath (pyresample.test.test_geometry.Test) ... ok
test_swath_equal (pyresample.test.test_geometry.Test) ... ok
test_swath_equal_area (pyresample.test.test_geometry.Test) ... ok
test_swath_not_equal (pyresample.test.test_geometry.Test) ... ok
test_swath_not_equal_area (pyresample.test.test_geometry.Test) ... ok
test_swath_wrap (pyresample.test.test_geometry.Test) ... ok
test_from_latlon (pyresample.test.test_grid.Test) ... ok
test_generate_linesample (pyresample.test.test_grid.Test) ... ok
test_latlons (pyresample.test.test_grid.Test) ... ok
test_latlons_mp (pyresample.test.test_grid.Test) ... ok
test_linesample (pyresample.test.test_grid.Test) ... ok
test_linesample_multi (pyresample.test.test_grid.Test) ... ok
test_proj4_string (pyresample.test.test_grid.Test) ... ok
test_proj_coords (pyresample.test.test_grid.Test) ... ok
test_resampled_image (pyresample.test.test_grid.Test) ... ok
test_resampled_image_mp (pyresample.test.test_grid.Test) ... ok
test_single_lonlat (pyresample.test.test_grid.Test) ... ok
test_image (pyresample.test.test_image.Test) ... ok
test_image_segments (pyresample.test.test_image.Test) ... /usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
ok
test_masked_image (pyresample.test.test_image.Test) ... ok
test_masked_image_fill (pyresample.test.test_image.Test) ... ok
test_nearest_neighbour (pyresample.test.test_image.Test) ... ok
test_nearest_neighbour_multi (pyresample.test.test_image.Test) ... ok
test_nearest_neighbour_multi_preproc (pyresample.test.test_image.Test) ... ok
test_nearest_resize (pyresample.test.test_image.Test) ... ok
test_nearest_swath (pyresample.test.test_image.Test) ... ok
test_nearest_swath_segments (pyresample.test.test_image.Test) ... /usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
ok
test_return_type (pyresample.test.test_image.Test) ... ok
test_area_con_cartesian_reduce (pyresample.test.test_kd_tree.Test) ... ok
test_area_con_reduce (pyresample.test.test_kd_tree.Test) ... ok
test_cartesian_reduce (pyresample.test.test_kd_tree.Test) ... ok
test_custom (pyresample.test.test_kd_tree.Test) ... ok
test_custom_base (pyresample.test.test_kd_tree.Test) ... ok
test_custom_multi (pyresample.test.test_kd_tree.Test) ... ok
test_custom_multi_from_sample (pyresample.test.test_kd_tree.Test) ... ok
test_custom_uncert (pyresample.test.test_kd_tree.Test) ... ok
test_gauss (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_base (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_fwhm (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_multi (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_multi_mp (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_multi_mp_segments (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_multi_mp_segments_empty (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_multi_uncert (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_sparse (pyresample.test.test_kd_tree.Test) ... ok
test_gauss_uncert (pyresample.test.test_kd_tree.Test) ... ok
test_masked_fill_float (pyresample.test.test_kd_tree.Test) ... ok
test_masked_fill_int (pyresample.test.test_kd_tree.Test) ... ok
test_masked_full (pyresample.test.test_kd_tree.Test) ... ok
test_masked_full_multi (pyresample.test.test_kd_tree.Test) ... ok
test_masked_gauss (pyresample.test.test_kd_tree.Test) ... ok
test_masked_multi_from_sample (pyresample.test.test_kd_tree.Test) ... ok
test_masked_nearest (pyresample.test.test_kd_tree.Test) ... ok
test_masked_nearest_1d (pyresample.test.test_kd_tree.Test) ... ok
test_nearest (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_1d (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_base (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_empty (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_empty_masked (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_empty_multi (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_empty_multi_masked (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_from_sample (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_mp (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_multi (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_multi_unraveled (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_remap (pyresample.test.test_kd_tree.Test) ... ok
test_nearest_segments (pyresample.test.test_kd_tree.Test) ... /usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1746: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  res = empty((N,)+dimensions, dtype=dtype)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1749: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  tmp.shape = (1,)*i + (dim,)+(1,)*(N-i-1)
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1751: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  val = zeros(newdim, dtype)
ok
test_reduce (pyresample.test.test_kd_tree.Test) ... ok
test_reduce_boundary (pyresample.test.test_kd_tree.Test) ... ok
test_area_def2basemap (pyresample.test.test_plot.Test) ... ok
test_easeplot (pyresample.test.test_plot.Test) ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
test_ellps2axis (pyresample.test.test_plot.Test) ... ok
test_orthoplot (pyresample.test.test_plot.Test) ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
test_plate_carreeplot (pyresample.test.test_plot.Test) ... ERROR
Testing if a point is inside an area. ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
Test how much two areas overlap. ... ok
Test if two areas overlap. ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
Testing the angle value between two arcs. ... ok
Test if two arcs intersect. ... ok
test_self_map (pyresample.test.test_swath.Test) ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
test_self_map_multi (pyresample.test.test_swath.Test) ... /home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py:76: UserWarning: All geometry objects expect longitudes in the [-180:+180[ range. We will now automatically wrap your longitudes into [-180:+180[, and continue. To avoid this warning next time, use routine utils.wrap_longitudes().
  'To avoid this warning next time, use routine utils.wrap_longitudes().')
ok
test_area_parser (pyresample.test.test_utils.Test) ... ok
test_load_area (pyresample.test.test_utils.Test) ... ok
test_not_found_exception (pyresample.test.test_utils.Test) ... ok
test_wrap_longitudes (pyresample.test.test_utils.Test) ... ok

======================================================================
ERROR: test_grid_filter (pyresample.test.test_geometry.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/test/test_geometry.py", line 444, in test_grid_filter
    (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py", line 514, in __init__
    proj = _spatial_mp.Proj(**proj_dict)
  File "/usr/lib/python2.7/dist-packages/pyproj/__init__.py", line 343, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1190)
RuntimeError: major axis or radius = 0 or not given

======================================================================
ERROR: test_grid_filter2D (pyresample.test.test_geometry.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/test/test_geometry.py", line 479, in test_grid_filter2D
    (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py", line 514, in __init__
    proj = _spatial_mp.Proj(**proj_dict)
  File "/usr/lib/python2.7/dist-packages/pyproj/__init__.py", line 343, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1190)
RuntimeError: major axis or radius = 0 or not given

======================================================================
ERROR: test_grid_filter_valid (pyresample.test.test_geometry.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/test/test_geometry.py", line 419, in test_grid_filter_valid
    (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py", line 514, in __init__
    proj = _spatial_mp.Proj(**proj_dict)
  File "/usr/lib/python2.7/dist-packages/pyproj/__init__.py", line 343, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1190)
RuntimeError: major axis or radius = 0 or not given

======================================================================
ERROR: test_latlong_area (pyresample.test.test_geometry.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/test/test_geometry.py", line 549, in test_latlong_area
    [-180, -90, 180, 90])
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py", line 514, in __init__
    proj = _spatial_mp.Proj(**proj_dict)
  File "/usr/lib/python2.7/dist-packages/pyproj/__init__.py", line 343, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1190)
RuntimeError: major axis or radius = 0 or not given

======================================================================
ERROR: test_plate_carreeplot (pyresample.test.test_plot.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/test/test_plot.py", line 53, in test_plate_carreeplot
    'test_files', 'areas.cfg'), 'pc_world')[0]
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/utils.py", line 109, in parse_area_file
    area_content)
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/utils.py", line 153, in _create_area
    config['AREA_EXTENT'])
  File "/home/antonio/debian/git/build-area/pyresample-1.1.3/pyresample/geometry.py", line 514, in __init__
    proj = _spatial_mp.Proj(**proj_dict)
  File "/usr/lib/python2.7/dist-packages/pyproj/__init__.py", line 343, in __new__
    return _proj.Proj.__new__(self, projstring)
  File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1190)
RuntimeError: major axis or radius = 0 or not given

----------------------------------------------------------------------
Ran 103 tests in 65.963s

FAILED (errors=5)

get_area_def balks if proj4_args is unicode

If the proj4 string passed to utils.get_area_def() is of type unicode instead of type str then the user gets RuntimeError: no projection control parameters specified.

Seems to me either _get_proj4_args() should try to handle the unicode or a more helpful error should be raised.

This is the string pulled from a netCDF file that made get_area_def() barf:
u'+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45'

Conda package and the use of the IOOS channel in the CI

Hi there, I am the maintainer of the IOOS channel and I noticed that pyresample is using that channel to get basemap on Windows. Note that the IOOS channel is now deprecated and our efforts are to maintain a better basemap recipe in the conda-forge channel. See https://github.com/conda-forge/basemap-feedstock and https://github.com/conda-forge/basemap-data-hires-feedstock

PS: I also created conda packages for pyresample here. Anyone interested in becoming a maintainer there is welcome to send a PR adding their GitHub handles to the maintainers list.

Change bilinear image test to use less memory

The bilinear resampling tests seem to use more memory than appveyor is comfortable with on the 32-bit Python 2.7 environment. We get occasional failed tests with the error:

======================================================================
ERROR: test_bilinear_multi (pyresample.test.test_image.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "c:\projects\pyresample-ly2q0\pyresample\test\test_image.py", line 249, in test_bilinear_multi
    area_con = msg_con.resample(self.area_def)
  File "c:\projects\pyresample-ly2q0\pyresample\image.py", line 402, in resample
    segments=self.segments)
  File "c:\projects\pyresample-ly2q0\pyresample\bilinear\__init__.py", line 85, in resample_bilinear
    epsilon=epsilon)
  File "c:\projects\pyresample-ly2q0\pyresample\bilinear\__init__.py", line 237, in get_bil_info
    _get_bounding_corners(in_x, in_y, out_x, out_y, neighbours, idx_ref)
  File "c:\projects\pyresample-ly2q0\pyresample\bilinear\__init__.py", line 426, in _get_bounding_corners
    y_diff = out_y_tile - in_y
MemoryError

It may be possible to change the test to not have large output or input data sizes.

Possibly inefficient array segmenting

I was looking through some of the nearest neighbor and kdtree stuff and found this section of code:

https://github.com/pytroll/pyresample/blob/master/pyresample/kd_tree.py#L337

    if segments > 1:
    # Iterate through segments
    for i, target_slice in enumerate(geometry._get_slice(segments,
                                                         target_geo_def.shape)):

        # Query on slice of target coordinates
        next_voi, next_ia, next_da = \
            _query_resample_kdtree(resample_kdtree, source_geo_def,
                                   target_geo_def,
                                   radius_of_influence, target_slice,
                                   neighbours=neighbours,
                                   epsilon=epsilon,
                                   reduce_data=reduce_data,
                                   nprocs=nprocs)

        # Build result iteratively
        if i == 0:
            # First iteration
            valid_output_index = next_voi
            index_array = next_ia
            distance_array = next_da
        else:
            valid_output_index = np.append(valid_output_index, next_voi)
            if neighbours > 1:
                index_array = np.row_stack((index_array, next_ia))
                distance_array = np.row_stack((distance_array, next_da))
            else:
                index_array = np.append(index_array, next_ia)
                distance_array = np.append(distance_array, next_da)

This code is saving memory by processing the data in chunks, but the end result will be the same size. So it is probably faster to allocate the entire end result arrays and write to the appropriate chunk rather than appending (which reallocates and copies the previous array).

Coordinates in geostationary projection

I tried to reproduce the geographical coordinates of an MSG/SEVIRI image using pyresample.geometry.AreaDefinition.get_lonlats() but found small differences from the results obtained by the CGMS formulas for the geostationary projection. Here is an example snippet:

import numpy as np
import matplotlib.pyplot as plt
import pyresample.geometry
import pyresample.plot

# Constants
RAD2DEG = 180.0/np.pi
DEG2RAD = np.pi/180.0
EQUATOR_EARTH_RADIUS = np.float64(6378.1370)
POLE_EARTH_RADIUS = np.float64(6356.7523)
SAT_ALTITUDE = np.float64(35785.863)  # above the equator
DIST_EARTHCENTRE_SAT = EQUATOR_EARTH_RADIUS + SAT_ALTITUDE  # = 42164
COFF = 1856
LOFF = 1856
NROWS = 3712
NCOLS = 3712
CFAC = -13642337*RAD2DEG
LFAC = -13642337*RAD2DEG
ANGULAR_PIXEL_RESOLUTION = np.float64(251.53) / 3.0  # in microrad


def image2scanangle(col, row):
    """
    Convert SEVIRI image coordinates to scanning angles (also called 
    'intermediate' coordinates in the CGMS document) corresponding to the
    centre of the image pixel.

    See appendix E in
    http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05057_SPE_MSG_LRIT_HRI&RevisionSelectionMethod=LatestReleased&Rendition=Web
    """
    x = (col - COFF) * ANGULAR_PIXEL_RESOLUTION/1E6
    y = (row - LOFF) * ANGULAR_PIXEL_RESOLUTION/1E6
    return x, y


def image2geo(col, row, lon_0=0.0):
    """
    Convert SEVIRI image coordinates to geographical lat/lon coordinates.

    See the "LRIT/HRIT Global Specification" document by CGMS, sections 4.4.3.2
    and 4.4.4:
    http://www.cgms-info.org/documents/cgms-lrit-hrit-global-specification-(v2-8-of-30-oct-2013).pdf

    @param col: SEVIRI image column(s) (Column 0 is West)
    @param row: SEVIRI image row(s) (Row 0 is North)
    @param lon_0: Nominal longitude (degrees east)
    @return: longitude(s), latitude(s)
    """
    # Compute scanning angles (or "intermediate coordinates")
    x, y = image2scanangle(col=col, row=row)

    # Transform scanning angles to geographical coordinates using the
    # geostationary projection (section 4.4.3.2)
    cosx = np.cos(x)
    cosy = np.cos(y)
    sinx = np.sin(x)
    siny = np.sin(y)

    # Some explanations to the numbers in the CGMS formulas:
    
    # a) The number 42164 in the formula for s_n is the distance between earth
    # centre and satellite
    
    # b) The number 1.006739501 in the formula for s_n and s_d is a 
    # "squared ellipsoid flatness" (EQUATOR_EARTH_RADIUS/POLE_EARTH_RADIUS)**2
    f2 = (EQUATOR_EARTH_RADIUS/POLE_EARTH_RADIUS)**2

    # c) I didn't find out where the number 1737122264 in the formula for 
    # s_d comes from...
    num = 1737122264

    sd = np.sqrt((DIST_EARTHCENTRE_SAT*cosx*cosy)**2 - (cosy**2 + f2*siny**2)*num)
    sn = (DIST_EARTHCENTRE_SAT*cosx*cosy - sd) / (cosy**2 + f2*siny**2)
    s1 = DIST_EARTHCENTRE_SAT - sn*cosx*cosy
    s2 = sn*sinx*cosy
    s3 = -sn*siny
    sxy = np.sqrt(s1**2 + s2**2)

    lon = RAD2DEG*np.arctan(s2/s1) + lon_0
    lat = RAD2DEG*np.arctan(f2*s3/sxy)

    return np.ma.masked_invalid(lon), np.ma.masked_invalid(lat)

if __name__ == '__main__':
    # Define image segment (full image here)
    col = np.arange(0, NCOLS)
    row = np.arange(0, NROWS)
    mcol, mrow = np.meshgrid(col, row)

    # Compute coordinates using CGMS formula
    lon, lat = image2geo(col=mcol, row=mrow, lon_0=0)

    # Compute coordinates using pyresample...
    #
    # ... determine projection corners. Relation between scanning angles
    # and projection coordinates is from: http://proj4.org/projections/geos.html
    x, y = image2scanangle(col=col, row=row)
    x_proj = x * SAT_ALTITUDE*1E3
    y_proj = y * SAT_ALTITUDE*1E3
    mycorners = [x_proj[0], -y_proj[-1], x_proj[-1], -y_proj[0]]  # The most eastern pixel is n=1856-1 and the most southern line is m=1856-1
    
    # ... setup AreaDefinition
    print mycorners
    seviri = pyresample.geometry.AreaDefinition(
        area_id='seviri',
        name='seviri',
        proj_id='seviri',
        proj_dict=dict(proj='geos',
                       lon_0=0,
                       h=SAT_ALTITUDE*1000,
                       a=EQUATOR_EARTH_RADIUS*1000,
                       b=POLE_EARTH_RADIUS*1000),
        x_size=NCOLS,
        y_size=NROWS,
        area_extent=mycorners
    )
    
    # ... compute coordinates and mask space pixels as they have invalid
    # values (1e30 or so)
    pyr_lon, pyr_lat = seviri.get_lonlats(dtype='float64')
    pyr_lon = np.ma.masked_where(np.logical_or(pyr_lon < -90, pyr_lon > 90),
                                 pyr_lon)
    pyr_lat = np.ma.masked_where(np.logical_or(pyr_lat < -90, pyr_lat > 90),
                                 pyr_lat)
    
    # Compute differences and plot results
    plt.figure()
    bmap = pyresample.plot.area_def2basemap(seviri)
    bmap.drawcoastlines()

    plt.figure()
    plt.imshow((pyr_lon - lon), cmap='bwr', vmin=-0.1, vmax=0.1)
    plt.title('lon diff')
    plt.colorbar()

    plt.figure()
    plt.imshow((pyr_lat - lat), cmap='bwr', vmin=-0.1, vmax=0.1)
    plt.title('lat diff')
    plt.colorbar()

    plt.show()

Did I understand/implement the CGMS formulas incorrectly or is that a pyresample issue? I'm using pyresample-1.5.0 .

commit a70cd80 broke scipy compatibility

@bhawkins, your pull request broke scipy compatibility: the data_pts attribute isn't available in the KDTree object.

DOC: explain how to trigger use of pykdtree ?

From [email protected] on May 31, 2013 12:02:43

Hei Esben,

We are running into some CPU-time bottleneck with pyresample (0.8.0).

We would like to try the new pykdtree support. But how do we do that?

Are there some installation flags to setup.py, for example?

I found nothing in the doc, but maybe I should look harder...

Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=11

multiprocessing installed with 'pip install' although not needed as in core of 2.7

From [email protected] on August 29, 2011 13:18:12

Hi,

this might as well be a problem with pip... What steps will reproduce the problem? 1. pip install pyresample-0.7.8.tar.gz What is the expected output? What do you see instead? The version of Python I use (2.7.1) also includes the multiprocessing module in its core distribution. I would expect that only missing modules are installed during the installation of pyresample. However, with pyresample-0.7.8, the "multiprocessing" backport for python 2.4 and 2.5 is fetched from pypi and installed. What version of the product are you using? On what operating system? Python 2.7.1, pip 0.8.1, distribute 0.6.14, OpenSUSE 11.1, python build with Intel compilers (so not the OpenSUSE python distribution) Please provide any additional information below. This might be a bug in pip (perhaps it doesn't check for core modules), but I thought I report it here first. IMHO, rather than installing it automatically, pyresample should only state the requirement and let the user decide to install a package which overwrites a core one. Or properly check the python version before making an attempt to install multiprocessing.

Original issue: http://code.google.com/p/pyresample/issues/detail?id=5

check usage of GridDefinition and longitude wrapping

From [email protected] on December 11, 2013 09:48:59

Check http://earthpy.org/interpolation_between_grids_with_pyresample.html 1. it was necessary to use SwathDefinition for orig_def and targ_def while I would have expected to use GridDefinition (using GridDefinition instead of SwathDefinition breaks the resampled field).

I had to substract 180 to the NCEP longitude (or the resampled was broken). That is not intuitive since exerything should happen in cartesian space.

Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=14

Unit test fail in test_swath for test_self_map and test_self_map_multi

I see two unit test fails from pyresample on Mac OS X 10.10 with numpy 1.9.1:

Python 2: https://gist.github.com/cdeil/f89d37a9314d518228a8
Python 3: https://gist.github.com/cdeil/c0767a434fafbe7a04d0

I think the fail is the same, but for some reason only on Python the values are printed by nose!???

======================================================================
FAIL: test_self_map (pyresample.test.test_swath.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/deil/code/pyresample/pyresample/test/test_swath.py", line 51, in test_self_map
    msg='Failed self mapping swath for 1 channel')
AssertionError: 668848.14481726696 != 668848.082208 within 1 places : Failed self mapping swath for 1 channel

======================================================================
FAIL: test_self_map_multi (pyresample.test.test_swath.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/deil/code/pyresample/pyresample/test/test_swath.py", line 71, in test_self_map_multi
    msg='Failed self mapping swath multi for channel 1')
AssertionError: 668848.14481726696 != 668848.082208 within 1 places : Failed self mapping swath multi for channel 1

----------------------------------------------------------------------

There's also a bunch of warnings that I guess should be caught in the library or tests if it's normal and expected behaviour.

Python 3.3 support?

From [email protected] on January 25, 2014 08:49:13

Are there plans to support python 3?

Original issue: http://code.google.com/p/pyresample/issues/detail?id=15

Can not install on Centos 7

I try to install on last Centos version and receive the following message:
"RuntimeError: maximum recursion depth exceeded in cmp"

See full output of install in the file attached
pyresamle_install_error.txt

Thank you.

Unittest fail on big-endian architectures

It seems that the test test_get_array_hashable fails on big endian architectures.
You can find all details about the conditions under which the test has been run at [1].

We are trying to include the new pyresample v1.7.1 in debian but this issue is blocking for us.
Can you pleas provide some hint about a possible solution/workaround?

[1] https://buildd.debian.org/status/package.php?p=pyresample

======================================================================
FAIL: test_get_array_hashable (pyresample.test.test_geometry.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/<<PKGBUILDDIR>>/pyresample/test/test_geometry.py", line 226, in test_get_array_hashable
    geometry.get_array_hashable(arr))
  File "/usr/lib/python2.7/dist-packages/numpy/testing/utils.py", line 1395, in assert_allclose
    verbose=verbose, header=header, equal_nan=equal_nan)
  File "/usr/lib/python2.7/dist-packages/numpy/testing/utils.py", line 778, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

(mismatch 50.0%)
 x: array([ 51,  51,  51,  51,  51,  51, 243,  63, 205, 204, 204, 204, 204,
       204, 244,  63, 102, 102, 102, 102, 102, 102, 246,  63,   0,   0,
         0,   0,   0,   0, 248,  63], dtype=uint8)
 y: array([ 63, 243,  51,  51,  51,  51,  51,  51,  63, 244, 204, 204, 204,
       204, 204, 205,  63, 246, 102, 102, 102, 102, 102, 102,  63, 248,
         0,   0,   0,   0,   0,   0], dtype=uint8)

----------------------------------------------------------------------
Ran 158 tests in 851.688s

FAILED (failures=1)

Negative Polygon Area

spherical_geometry.get_polygon_area() computes negative values near the poles. The result also depends on the order of the given corners:

from pyresample.spherical_geometry import Coordinate, get_polygon_area

ll = Coordinate(0, 89.925)
ul = Coordinate(0, 89.975)
ur = Coordinate(0.05, 89.975)
lr = Coordinate(0.05, 89.925)


print get_polygon_area([ll, ul, ur, lr])  # -0.00043633152367
print get_polygon_area([ll, ul, lr, ur])  # 9.96988713808e-10

Is there an expected order of the corners?

Packaging a .egg-info with pyresample?

From [email protected] on May 06, 2013 08:33:05

Hi,

This is more a question than a defect. When I "python setup.py install" other modules than pyresample, it is not rare to obtain a ".egg-info" file installed under site-packages/. This is not the case for pyresample:

Here is an example from a software we have been developing and that builds several packages from source:
<kapphorn|pmw_decode> =:D ls mypython/lib/python2.6/site-packages/*.egg-info
mypython/lib/python2.6/site-packages/basemap-1.0-py2.6.egg-info
mypython/lib/python2.6/site-packages/bufr_metadb-0.1-py2.6.egg-info
mypython/lib/python2.6/site-packages/bufr_transform-0.1-py2.6.egg-info
mypython/lib/python2.6/site-packages/lockfile-0.9.1-py2.6.egg-info
mypython/lib/python2.6/site-packages/netCDF4-0.9.1-py2.6.egg-info
mypython/lib/python2.6/site-packages/python_bufr-0.2_3-py2.6.egg-info
mypython/lib/python2.6/site-packages/setuptools-0.6c11-py2.6.egg-info
mypython/lib/python2.6/site-packages/sqlalchemy_marshall-1.0-py2.6.egg-info
mypython/lib/python2.6/site-packages/xml_marshall-1.0-py2.6.egg-info

BUT:
ls mypython/lib/python2.6/site-packages/pyresample-0.7.7-py2.6.egg-info
no such file or directory.

So I wonder if this is common practice to have this .egg-info installed and why does not appear for pyresample?

I use the .egg-info as a target for my Makefile rules: if the .egg-info is newer than the tar-ball containing the source code, do not "python setup.py install". And I cannot apply the same trick for pyresample.

Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=9

Join user and package area files

Currently the parse_area_file function only allows one single file. This causes limitations in software like satpy where we may want to have a "user" area file and a "package" area file. The easiest/cleanest way to handle this is to make a pyresample allow more than one file. This is even more useful for the new yaml files which could allow users area files to "update" existing information in the "package" area file.

I propose that there be a new function (for backwards compatibility) that accepts one of the following input formats for the areas files:

single area file name
single file-like object
single file content string
iterable of area file name or file-like object

Correct resampling of multiple swaths that do not fit into memory

We are currently working on the problem of resampling multiple swaths in a correct manner. Our main concern are edge effects when stitching swaths or parts of orbits together. We mostly use a hamming window for resampling so we do not want to ignore the edges of our swaths. We do have a working solution for the problem and I am currently looking into the best way to open source it.

My main questions are:

Would you consider this functionality a part of pyresample or should wrappers around pyresample manage inter-swath resampling.
Are you aware of other projects that tackle the same issue which we could contact if your answer to question 1 is no.

Otherwise we will open source it as our own project and keep you informed if you are interested.

Computing density_of_x (alternatively "counting number of x)" while re-gridding

From [email protected] on May 06, 2013 08:42:38

Hi,

This is more a feature request:

I have a fine resolution (1/20 deg) global land/water mask. Each grid cell has thus a lat, lon, and a byte code for "ocean", "land", "lake", etc...

I would like to compute the density_of_ocean, density_of_land, density_of_lake, etc... on a, say, polar EASE grid with 25 km spacing. A good approximation for my density_of_land is the ratio between the number of "land" original grid cells falling into a given "target grid" cell, to the total number of original grid cells falling into that same "target grid" cell.

Does pyresample have such a feature?

An alternative for me would be to first remap (nearest neighbour) the 1/20 deg map to a, say 2.5km EASE grid, with grid cells perfectly aligned with my target (25km) grid. Then I can do the density/counting myself.

But I am interested to hear if pyresample can/could do that out of the box :)

Cheers,
Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=10

release 1.1.1 breaks API

I'm not sure how pyresample should be imported but i used it as shown in the example below up until version 1.1.0.

For example:

>> import pyresample as pr
>>> pr.__version__
'1.1.1'
>>> pr.kd_tree.which_kdtree()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'module' object has no attribute 'kd_tree'
>>> pr.geometry.SwathDefinition
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'module' object has no attribute 'geometry'

Test failures on mips*, hppa & hurd-i386

We update the pyresample debian package to v1.2.4 and now we get build failures on mips*, hppa & hurd-i386.

The error looks like:

FAIL: test_gauss_uncert (pyresample.test.test_kd_tree.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/«PKGBUILDDIR»/pyresample/test/test_kd_tree.py", line 117, in test_gauss_uncert
    len(w) != 1, 'Failed to create neighbour warning')
AssertionError: Failed to create neighbour warning

----------------------------------------------------------------------
Ran 114 tests in 582.981s

FAILED (failures=3)
E: pybuild pybuild:274: test: plugin distutils failed with: exit code=1: python2.7 setup.py test 
dh_auto_test: pybuild --test -i python{version} -p 2.7 returned exit code 13
debian/rules:14: recipe for target 'build-arch' failed
make: *** [build-arch] Error 25

Complete log available at

https://buildd.debian.org/status/fetch.php?pkg=pyresample&arch=mips&ver=1.2.4-2&stamp=1468065945
https://buildd.debian.org/status/fetch.php?pkg=pyresample&arch=mipsel&ver=1.2.4-2&stamp=1468066763
https://buildd.debian.org/status/fetch.php?pkg=pyresample&arch=hppa&ver=1.2.4-2&stamp=1468067276
https://buildd.debian.org/status/fetch.php?pkg=pyresample&arch=hurd-i386&ver=1.2.4-2&stamp=1468065283
https://buildd.debian.org/status/fetch.php?pkg=pyresample&arch=mips64el&ver=1.2.4-2&stamp=1468065703

Can you please provide some hint on how to fix it?

Using GridDefinition with GOES-12 data

I'm trying to get an area definition using GridDefinition for GOES-12 data. But as you can see in the image below, the lats and lons in the top left corner are not on the sphere of the Earth so I'm getting the error:

Some latitudes are outside the [-90.;+90.] validity range.

Is there a work around for this or will I have to use AreaDefinition?

Preprocces Area Def to grid Def

I'd like to be able to construct an AreaDefinition and GridDefinition and preprocess the grid resampling like so

https://pyresample.readthedocs.io/en/latest/preproc.html

but to have one be a Grid definition instead of a AreaDefinition.

Is this possible?

Test suite stalls

The test suite stalls on debian sid amd64 with python3.5:

$ pybuild base:184: python3.5 setup.py test 

...

test_nearest_resize (pyresample.test.test_image.Test) ...

Please note that the test seems to pass regularly with python 2.7.

Patch for importing version in /trunk/setup.py

From jakob.malm on August 04, 2011 11:48:31

Import version directly using imp.load_source. This way, any previous imports of a module named "version" will not affect this import. In addition, the sys.path variable remains unchanged. And, as before, the pyresample package is not imported.

Attachment: setup.py.patch

Original issue: http://code.google.com/p/pyresample/issues/detail?id=4

Patch for /setup.py

From jakob.malm on November 01, 2011 18:52:57

Check python version before adding multiprocessing to list of required packages. Otherwise, the package gets downloaded and installed despite it being in the standard library for Python >= 2.6.

(Another way to do it could of course be to try to import multiprocessing. I don't know which would be preferrable.)

Attachment: setup.py.patch

Original issue: http://code.google.com/p/pyresample/issues/detail?id=6

Missing basemap and matplotlib dependencies

Basemap and matplotlib are required for pyresample/plot but are not specified as a requirement or extra, so setuptools don't know about them and test_plot.py fails in a local vanilla environment.

Proposal: New area file format

This is a proposal for pyresample to support a new file format to get area definitions from. There are a lot of decisions that could go in to something like this so I wanted to make an issue before I threw something together and have everyone hate it. Depending on what gets decided the current format may stay around, but my hope is that the new format is so good no one will want to use the old format anymore. Also this new format may be made most flexible if there is a method or function or class that provides most of the convenience that the configuration provides (see below). What this would mean for the current classes and methods, I'm not sure.

I've been thinking about this for a while now during satpy development because satpy will be heavily used by Polar2Grid and therefore Polar2Grid's users. Polar2Grid has its own comma-separated way of specifying "grids" and I'd like to use whatever format satpy uses, but I'm not really a fan of the current format nor does the current format provide all of the features that P2G's provides. These features revolve around the idea of convenience for the user and allowing them to enter what they want or what they know and having logical defaults fill in the gaps. In P2G's case, some users don't know what projections are or don't know that meters in LCC are different than meters in Mercator.

NOTE: All areas have a name and an optional description. The name must be unique among all other areas in the file.

Ways Areas could be specified but not limited just to these:

Extents - Outer pixel edge coordinates in meters (legacy spec)
a. PROJ.4 String (or dictionary, debatable)
b. 4-element list of extents in projection units (same as current format)
c. number of columns and rows in pixels
Geotiff specification
a. PROJ.4 string or dict
b. upper-left origin (X, Y) in projection units
c. Cell size (X, Y) in projection units (a.k.a. pixel resolution)
d. number of columns and rows in pixels
Circle - resulting area is still rectangular
a. PROJ.4 string or dict
b. center coordinate (X, Y) in projection units
c. radius in projection units
Area of interest - Geotiff based on center point
a. Center coordinate (X, Y) in projection units
b. Cell size (X, Y) in projection units
c. number of columns and rows in pixels
d. PROJ.4 string or dict
Swath Definition - Maybe?
a. binary longitude file (.npz is numpy array, .dat is float32 binary, etc?)
b. binary latitude file

In addition to this I think it should be possible for a user to specify anything that is in "projection units" in degrees except for maybe cell size and radius. I've often helped people who started using pyresample and one of their first questions is "What are these 4 numbers for extents? How do I calculate them?". Allowing for someone to enter degrees would help with that.

Another feature would be logical defaults. For example, in number 3 or 4 the PROJ.4 could be optional and a "best projection" could be chosen based on the absolute value of the latitude of the center point (ex. mercator for < 15, lcc for <60, polar stereographic otherwise).

Another feature that would require support inside of pyresample would be "dynamic grids" as they are called in Polar2Grid. These are grids where things like origin or number of cols/rows or even cell size is determined by the data being resampled. Based on my experience in Polar2Grid with this an important thing that can be overlooked is that datasets that are being resampled for the same geographic scene (time and location) should usually share the same resolved dynamic Area. Most likely this resolved Area used the highest resolution input swath definition to determine things. This is very important if the resampled results are going to be compared or merged in to a composite. In Polar2Grid this resolution step of dynamic areas is hidden inside the resampling code. This probably wouldn't work best for pyresample since it limits what the user can do with the resulting area.

Lastly, what format would best serve something like this:

YAML (my vote)
INI
Other?

python 3 compatibility broken with unicode patch

Support for lat-long swath data

From [email protected] on October 12, 2012 10:13:43

It would be nice to include support for swath data in lat-long (EPSG:4326). Currently, if you feed in area_extent in degrees and set projection information to EPSG:4326, it doesn't quite work.

Maybe using is_latlong() flag from the result from proj4 can help resolve this?

Piyush

Original issue: http://code.google.com/p/pyresample/issues/detail?id=7

utils.parse_area_file does not work as expected

From [email protected] on May 18, 2011 11:45:20

Hei,

I am trying to use utils.parse_area_file('file',['area_1', 'area_2']) as described in the function documentation (from ipython):

Type: function
Base Class: <type 'function'>
String Form: <function parse_area_file at 0x95aa224>
Namespace: Interactive
File: /home/thomasl/software/pmw_decode/mypython/lib/python2.6/site-packages/pyresample-0.7.6-py2.6.egg/pyresample/utils.py
Definition: pyresample.utils.parse_area_file(area_file_name, *regions)
Docstring:
Parse area information from area file

:Parameters:
area_file : str
    Path to area definition file
regions : str argument list 
    Regions to parse. If no regions are specified all 
    regions in the file are returned

:Returns:
area_defs : list
    List of AreaDefinition objects

So it says here that the second parameter can be a list of regions, and that the return value should then be a list of area_defs.

But:

res = utils.parse_area_file('file','area_1')
OK (res is a list object)

res = utils.parse_area_file('file',['area_1'])
AreaNotFound: Area "['area_1']" not found in file "../par/grids.def"

Cheers,
Thomas

Original issue: http://code.google.com/p/pyresample/issues/detail?id=2

Assigning to slices/views of MaskedArray generates a warning

The KD tree code manipulates the mask, which raises a numpy warning.

To replicate:

from satpy.scene import Scene
global_scene = Scene(platform_name = "SNPP", sensor = "viirs", 
    base_dir = "viirs_test_download_2017-08-25_ND")
global_scene.load(["M10"])
global_scene.resample("worldeqc3km70")

I'm using numpy 1.13.1, satpy 0.7.1 and pyresample 1.5.0.

Need the outer boundary of the area extent in lon,lat for all four corners

From [email protected] on February 26, 2014 13:48:36

Sometimes one needs to get the lon,lat of the foure corners of the area extent. Not the lon,lat of the four corner pixels, but the actual corners of the area.

I therefore added a new property to the AreaDefinition class, called "outer_boundary_corners". I did this in the pre-master branch

Usage:

In [19]: area.corners
Out[19]:
[(3.6330765317047953, 65.814387639583856),
(33.719973797301968, 64.899876035407289),
(27.673973672109508, 52.700325218767617),
(6.893594960144088, 53.335796768826384)]

In [20]: area.outer_boundary_corners
Out[20]:
[(3.6206470940779134, 65.818140630566901),
(33.733114510263896, 64.902858589601166),
(27.67965639582464, 52.694968994188152),
(6.8869748038970524, 53.330819380313962)]

I added a simple test also (for both corners and outer_boundary_corners).

Original issue: http://code.google.com/p/pyresample/issues/detail?id=16

TRMM Data Sample

Hello,

    I want to sample TRMM 2A25 correctZFactor data over passes from radar site.The radar takes  observation in 150 Km radius (72.7 <= longitude < = 74.7 and  17.5,<=latitude <= 19.5).

Longitude (nscan X rays)(9249 X 49)
Latitude (9249 X 49)
correctZFactor (9247 x 49 x 80 ((nscan X rays X height)
please help me to sample data

duplicate an AreaDefinition object with coarser/finer definition

Add a feature for the AreaDefinition class (https://pyresample.readthedocs.org/en/latest/geo_def.html#areadefinition). I would need a class method that creates a new AreaDef object from an existing one. The new AreaDef has the same projection and coordinate extent, but the grid_spacing (equivalently the number of gridpoints x_size and y_size) is scaled by a factor provided by the user. This would return a new AreaDef that matches exactly the existing one, but with twice (thrice, four times as much) the original resolution.

Example:

area_def = geometry.AreaDefinition(area_id, name, proj_id, proj_dict, x_size,y_size, area_extent)
new_area_def = area_def.finer_grid_spacing(3,2)

new_aera_def now has three times as much cells in the x axis, and twice in the y axis. If only 1 value is given, the same is used for x and y.

Extend custom weight functions

For some use cases the radius is not enough to do a proper resampling.

We are having this problem with categorical data which can not be averaged in a meaningful way but needs e.g. a decision tree to determine the correct aggregate.
This is probably slightly related to issue #12.

Area definition have empty projection_{x, y}_coords by default

From [email protected] on November 26, 2013 22:20:18

What steps will reproduce the problem? 1. create an area definition
2. area.projection_x_coords is None
3. calling area.get_proj_coords(cache=True) populates the attributes What is the expected output? What do you see instead? I would expect the attributes to be populated, at least when accessed. They are None now. What version of the product are you using? On what operating system? v1.1.0 Please provide any additional information below. calling area.get_proj_coords(cache=True) populates the attributes. Shouldn't it be done on the fly ?

Original issue: http://code.google.com/p/pyresample/issues/detail?id=13

resample area EPGS 32735

From [email protected] on November 10, 2012 21:39:05

I have defined an area in South Africa (WGS 84 UTM 35S) using the following args:
REGION: Douglas_SA_300m {
NAME: Douglas_SA_300m
PCS_ID: Douglas_SA_300m
PCS_DEF: proj=utm, zone=35, south, ellps=WGS84, datum=WGS84, units=m
XSIZE: 430
YSIZE: 280
AREA_EXTENT: (164000.0, 6750000.0, 293000.0, 6834000.0)
};

The proj4_args can´t load. If I delete (south)in the defined area the arguments load but the file is resampled to EPSG:32635 (WGS 84 UTM 35N)instead of 35S.

I would liek to know how to introduce the south zone in UTM on the area definition?

Ernesto

Original issue: http://code.google.com/p/pyresample/issues/detail?id=8

reduce_data bug with polar grids

Hello,

This is a bug I have been chasing for a while, but I have the impression it somehow is worse now (1.5.0) than before. It happens when using kd_tree.resample_* and reduce_data=True (which is the default) with polar orbiting satellites, on polar grids. At least the polar grids I am using for my work.

I prepared a script to generate test data and compare the remapping with and without reduce_data. It results in the image below.

To me, the reduce_data feature should only speed-up the resampling, not change the result in any way. So this is a serious flaw.

I hope this can be made rather easily into a unit test for reduce_data since it is self contained (no need for external data files to be loaded).

Where should we go from there? The github tool will not let me attach a .py. I can send it by mail. I try to paste it here but it does not look good.

import numpy as np
import pyresample as pr
import matplotlib.pylab as plt
from mpl_toolkits.basemap import Basemap
import pyproj

print "Pyresample {}".format(pr.__version__,)

# routine to create a polar area_def (NSIDC's 25km SSM/I grid)
def get_areadef(hemis):
    nx, ny = {'nh':(448,304), 'sh':(332,316)}[hemis][::-1]
    a = 25000
    if hemis == 'nh':
        xl = -3850000.00
        xL = xl+0.5*a + (nx-1)*a + 0.5*a
        yL = 5850000.00
        yl = yL-0.5*a - (ny-1)*a - 0.5*a
        pdict = {u'a': u'6378273', u'b': u'6356889.44891', u'lat_ts': u'70', u'lon_0': u'-45', u'proj': u'stere', u'lat_0': u'90', u'units':u'm'}
    elif hemis == 'sh':
        xl = -3950000.00
        xL = xl+0.5*a + (nx-1)*a + 0.5*a
        yL = +4350000.00
        yl = yL-0.5*a - (ny-1)*a - 0.5*a
        pdict = {u'a': u'6378273', u'b': u'6356889.44891', u'lat_ts': u'-70', u'lon_0': u'0', u'proj': u'stere', u'lat_0': u'-90', u'units':u'm'}
    else:
        raise ValueError("Not a valid area ('nh' and 'sh' are valid)")

    area_def  = pr.geometry.AreaDefinition('nsidc_ps_{}'.format(hemis,),
                                          'nsidc_ps_{}'.format(hemis,),
                                          'nsidc_ps_{}'.format(hemis,),
                                          pdict,
                                          nx, ny,
                                          [xl, yl, xL, yL])
    if area_def.pixel_size_x != 25000 or area_def.pixel_size_y != 25000:
        raise ValueError("The grid parameters do not match a grid spacing of 25km! (x:{} y:{})".format(area_def.pixel_size_x, area_def.pixel_size_y))
    
    return area_def


area_def = get_areadef('nh')

#simulate swath data:
#  1) have a single line of data going through (North) pole
slon = 125.
swath_lats_1 = np.arange(30.,90.,1.)
swath_lons_1 = np.array([slon,]*swath_lats_1.size)
swath_lats_2 = swath_lats_1[::-1]
swath_lons_2 = np.array([180+slon,]*swath_lats_2.size)

swath_lats = np.concatenate((swath_lats_1,[90.,],swath_lats_2))
swath_lons = np.concatenate((swath_lons_1,[0.,],swath_lons_2))

#  2) use pyproj to obtain x,y coordinates in the map
P = pyproj.Proj(area_def.proj4_string)
x,y = P(swath_lons, swath_lats)

#  3) use pyproj(inverse=True) to get 5 neighbouring lines that are side by side and not going through Pole. This is our simulated polar orbiting swath
nbscanpos = 5
lons = np.empty((swath_lons.size,nbscanpos))
lats = np.empty((swath_lats.size,nbscanpos))
data = np.empty((swath_lats.size,nbscanpos))
off = 100*1000  # be at least 100km from Pole
for i in range(0,nbscanpos):
    xp = x + off
    yp = y - off
    lons[:,i], lats[:,i] = P(xp,yp,inverse=True)
    data[:,i] = i+5
    off += 50*1000 # distance between scan positions is 50km


# Do 2 kdtree.resample() calls: exactly the same except for reduce_data=False|True (default is True)
swath_def = pr.geometry.SwathDefinition(lons, lats)
res_rdF = pr.kd_tree.resample_nearest(swath_def,data,area_def,radius_of_influence=50000,reduce_data=False)
res_rdT = pr.kd_tree.resample_nearest(swath_def,data,area_def,radius_of_influence=50000,reduce_data=True)

#plot
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(121)
ax.set_title('kdtree data_reduce=False')
bmap = pr.plot.area_def2basemap(area_def)
bmap.drawcoastlines()
bmap.imshow(res_rdF,interpolation='none',cmap=plt.cm.plasma, origin='upper')

ax = fig.add_subplot(122)
ax.set_title('kdtree data_reduce=True')
bmap = pr.plot.area_def2basemap(area_def)
bmap.drawcoastlines()
bmap.imshow(res_rdT,interpolation='none',cmap=plt.cm.plasma, origin='upper')


plt.savefig('./pyresample_reduce_data_bug.png',bbox_inches='tight')
print "./pyresample_reduce_data_bug.png ready"

plt.show()

pytroll / pyresample Goto Github PK

pyresample's Introduction

pyresample's People

Stargazers

Watchers

Forkers

pyresample's Issues

Recommend Projects

Recommend Topics

Recommend Org