tuw-geo / pygeogrids Goto Github PK

Grid shape is expected to be given as (size_lat, size_lon) when creating a grid:

When storing the grid, the shape is passed to the netcdf writer:

But later in save_lonlat() where it is used to structure the variables, latsize has to be the second, resp. lonsize has to be the first element, therefore it won't work.

Suggested fix: Change

latsize = global_attrs['shape'][0]
lonsize = global_attrs['shape'][1]

A numpy number is a single number with a numpy type

>>> a = np.arange(2)
>>> type(a[0])
<type 'numpy.int64'>

If indexed they throw a IndexError and not a TypeError like implemented here

Loading a grid from netcdf allows selecting a subset_flag to reduce the grid based on a single maskin the netcdf file. If we allow not only flags, i.e. subset==1, but variables, we could reduce grids from files directly based on e.g. landcover classes or climate classes (e.g. from esa cci landcover). Should be easy to implement, any objections?

IndexError should also be included in except statement to cover the possibility that gpi is given as numpy integer.

https://github.com/TUW-GEO/pygeogrids/blob/master/pygeogrids/grids.py#L663

This method is needed since type checking is error prone

it is not available in windows at the moment

set shape automatically from them

was created by pyscaffold and is not needed in out case

Retrieving lonlat information from a gpi fails in the function gpi2lonlat in this line , not sure if this is also related to #39

e.g. gp = 2412975 of the WARP grid cannot be found in if the activearrlon and activearrlat are just the land points.

If I set self.gpidirect to True, the function works.

It would be nice if one could specify the number of nearest neighbors to return when querying the kdtree. https://github.com/TUW-GEO/pygeogrids/blob/master/pygeogrids/nearest_neighbor.py#L191

The problem is that get_bbox_grid_points does not look a cell boundaries so it can not be used in the grid_points iterator directly.

Right now the output of
gpi2lonlat and gpi2rowcol do not fit together well since the row is the latitude and the column is the longitude.

in find_nearest_gpi

• grid generation,
• find nearest neighbor
• lookuptable generation
• storing and loading as netcdf file

I think we only build for python 2.7? I just tried to setup the develop env locally and using the environment.yml (same as in the CI) always install python2. See here: https://travis-ci.org/TUW-GEO/pygeogrids/jobs/546953517#L675

For some reason Basic Grid can not be pkled, but CellGrid can. This can lead to issues with eg parallelization frameworks and should be fixed.

import cloudpickle
from pygeogrids.grids import CellGrid, BasicGrid

cloudpickle.dumps(CellGrid([1,2,3], [1,2,3], [1,2,3]))
print("Cellgrid could be pickled")
cloudpickle.dumps(BasicGrid([1,2,3], [1,2,3]))

Traceback (most recent call last):
  File "/home/wpreimes/miniforge3/envs/smecv_gapfill/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3553, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-3-7762823bd9a7>", line 6, in <module>
    cloudpickle.dumps(BasicGrid([1,2,3], [1,2,3]))
  File "/home/wpreimes/miniforge3/envs/smecv_gapfill/lib/python3.10/site-packages/cloudpickle/cloudpickle.py", line 1479, in dumps
    cp.dump(obj)
  File "/home/wpreimes/miniforge3/envs/smecv_gapfill/lib/python3.10/site-packages/cloudpickle/cloudpickle.py", line 1245, in dump
    return super().dump(obj)
  File "<stringsource>", line 2, in pykdtree.kdtree.KDTree.__reduce_cython__
TypeError: no default __reduce__ due to non-trivial __cinit__
``

Otherwise stored and loaded grids are not equal

I was going to write my own grid management python package and I discovered this one, do you happen to publish the documentation somewhere?

I would be glad to contribute if there are synergies.

This module has been end-of-life'd and removed in favor of cartopy.

Depending on whether the NN search returns a distance containing a np.inf or not decides whether the returned indices and distances are a 1d array (mask applied) or a 2d array (no mask applied, i.e. no np.inf entry in the distances)

https://github.com/TUW-GEO/pygeogrids/blob/master/src/pygeogrids/grids.py#L468

Moved from TUW-GEO/pytesmo#9

other.kdTree.find_nearest_index in calc_lut returning no result with max_dist defined, results are returned without defining max_dist.

• finding nearest gpi
• getting cell numbers

Maybe define a API for slicing a grid

grid.find_nearest_gpi() takes a keyword max_dist, but this is ignored in the next line when calling ´´find_k_nearest_gpi()´´ and set to np.Inf

I guess L382 should be changed to max_dist=max_dist

Also the unit of max_dist is not clear from the docstring

save_lonlat does not handle arrlon/arrlat/gpi correctly if the order is different than the one that is used by numpy's reshape method. this problem also arises when especially the order of the gpi is different than the one that is presumed in this function, e.g. gpi is counting from "upper left" to "lower right".

load_grid with newer netcdf versions creates masked arrays, e.g. when calling grid.get_grid_points().
I think this is not intended and due to changes in the netcdf4 library in the back. We don't utilize masked arrays and in most cases they should be compatible, still should we fix that to restore the original behavior?

I get an error when using the new find_nearest_gpi function (with max_dist) from #66

>> import pygeogrids
>>pygeogrids.__version__
'0.3.1'
>> grid = pygeogrids.BasicGrid(lon=[16,17], lat=[45,46], gpis=[100,200])
>> grid.find_nearest_gpi(0,0)
(100, 5082354.193805456)
>> grid.find_nearest_gpi(0,0, max_dist=1000)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Users\wpreimes\Anaconda3\envs\ismn\lib\site-packages\pygeogrids\grids.py", line 400, in find_nearest_gpi
    gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=max_dist, k=1)
  File "C:\Users\wpreimes\Anaconda3\envs\ismn\lib\site-packages\pygeogrids\grids.py", line 441, in find_k_nearest_gpi
    gpi = self.activegpis[ind]
IndexError: arrays used as indices must be of integer (or boolean) type

I think the function should return an empty array instead.

It works when I dont set the gpis manually when creating the grid

>> grid = pygeogrids.BasicGrid(lon=[16,17], lat=[45,46])
>> grid.find_nearest_gpi(0,0, max_dist=1000)
(array([], dtype=float64), array([], dtype=float64))

At some point the GPI is used for looking up the lon/lat I think, which only works for default GPIs (starting at 0)

@s-scherrer can you fix this?

the methods gpi2cell gpi2lonlatusenumpy.searchsorted` but this function does not really search for the occurrence of the gpi in the array but just for the position of where to insert the gpi. This is fine as long as the gpi is actually in the grid. We need a better check/search algorithm for this.