If it takes another couple of years to break, and need another re-write that's probably OK.

But if there is a better/more standard way to do what I'm trying to do, please let me know.

The underlying issue is that I have a netcdf with a RotatedPole grid (e.g. any of the CORDEX data), and I need to find the nearest grid point to a certain point. To do this, I map my point onto the RotatedPole projection (converting it from lat, lon, to x, y, or rlat, rlon), and then use the index mapping's nearest method to find the nearest grid point.

It didn't seem like there was an obvious way to do this with the actual API functions, so that's what I landed on.

Thank you, this is sufficient for me for now. I was able to replace

python nearest_point = remap_label_indexers(self.data, dict(x=x, y=y), method='nearest')[0]

with

python nearest_point = map_index_queries(self.data, dict(x=x, y=y), method='nearest').dim_indexers

@shoyer That would either not work, or be needlessly expensive, I think. The message generation might be expensive (e.g. if I want a sum or mean of the differences). With a call back it only happens if it is needed. With a pre-computed message it would be computed every time.. Correct me if I'm wrong.

@TomNicholas My particular use case is that have datasets that are large enough that I can't see the full diff, so might miss major changes. I'm wanting to pass in something like lambda a, b: f"Largest difference in data is {abs(a-b).max().item()}", so I can quickly see if the changes are meaningful. Obviously a more complex function might also be useful, like a summary/describe table output of the differences..

I know I could set the tolerances higher, but the changes are not numerical errors, and I want to see them before updating the test data that they are comparing against.

Entirely possible that there are better ways to do this, of course :)

I'm seeing a problem with geotiffs that I think might be related:

sh $ gdalinfo /home/nedcr/cr/software/ana/lib/geo/tests/data/Urandangi_MGA55.tif Driver: GTiff/GeoTIFF Files: /home/nedcr/cr/software/ana/lib/geo/tests/data/Urandangi_MGA55.tif /home/nedcr/cr/software/ana/lib/geo/tests/data/Urandangi_MGA55.tif.aux.xml Size is 10, 10 Coordinate System is: PROJCS["GDA94 / MGA zone 55", GEOGCS["GDA94", DATUM["Geocentric_Datum_of_Australia_1994", SPHEROID["GRS 1980",6378137,298.257222101, AUTHORITY["EPSG","7019"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6283"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4283"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",147], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",10000000], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","28355"]] Origin = (-406507.954209543997422,7588834.152862589806318) Pixel Size = (263.500000000000000,-265.500000000000000) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( -406507.954, 7588834.153) (138d16' 6.99"E, 21d34'24.95"S) Lower Left ( -406507.954, 7586179.153) (138d16' 1.84"E, 21d35'50.30"S) Upper Right ( -403872.954, 7588834.153) (138d17'37.56"E, 21d34'29.73"S) Lower Right ( -403872.954, 7586179.153) (138d17'32.42"E, 21d35'55.09"S) Center ( -405190.454, 7587506.653) (138d16'49.70"E, 21d35'10.02"S) Band 1 Block=10x10 Type=Float32, ColorInterp=Gray Min=0.069 Max=8.066 Minimum=0.069, Maximum=8.066, Mean=2.556, StdDev=1.749 NoData Value=-3.40282306073709653e+38 Metadata: STATISTICS_MAXIMUM=8.06591796875 STATISTICS_MEAN=2.5563781395387 STATISTICS_MINIMUM=0.068740844726562 STATISTICS_STDDEV=1.7493082797107 STATISTICS_VALID_PERCENT=89

```python import xarray as xr from osgeo.gdal import Open

ras = Open('/home/nedcr/cr/software/ana/lib/geo/tests/data/Urandangi_MGA55.tif') ras.GetGeoTransform()

(-406507.954209544, 263.5, 0.0, 7588834.15286259, 0.0, -265.5)

ds = xr.open_rasterio('/home/nedcr/cr/software/ana/lib/geo/tests/data/Urandangi_MGA55.tif') ds.transform

(263.5, 0.0, -406507.954209544, 0.0, -265.5, 7588834.15286259)

```

The transform in the xarray dataset is transposed, and not really useful anymore.

Output of xr.show_versions()

I guess ds = ds.drop([c for c in ds.coords if c not in ds.dims]) works. Might be nice to have a convenience function for that though.

Is there currently a way to drop unused coordinates, like this?

Awesome, thanks shoyer! :)

Ok, no worries. I'll try it if it gets desperate :)

Thanks for your help, shoyer!

So it would be more efficient to concat all of the datasets (subset for the relevant variables), and then just use a single .to_dataframe() call on the entire dataset? If so, that would require quite a bit of refactoring on my part, but it could be worth it.

Yeah, I'm loading each file separately with xr.open_dataset(), since it's not really a multi-file dataset (it's a lot of single-site datasets, some of which have different variables, and overlapping time dimensions). I don't think I can avoid loading them separately...

Not easily - most scripts require multiple (up to 200, of which the linked one is one of the smallest, some are up to 10Mb) of these datasets in a specific directory structure, and rely on a couple of private python modules. I was just asking because I thought I might have been missing something obvious, but now I guess that isn't the case. Probably not worth spending too much time on this - if it starts becoming a real problem for me, I will try to generate something self-contained that shows the problem. Until then, maybe it's best to assume that xarray/pandas are doing the best they can given the requirements, and close this for now.

I loaded it from a netcdf file. There's an example you can play with at https://dl.dropboxusercontent.com/u/50684199/MitraEFluxnet.1.4_flux.nc

Slightly slower (using %timeit in ipython)

Squeeze is pretty much identical in efficiency. Seems very slightly better (2-5%) on smaller datasets. (I still need to add the final [data_vars] to get rid of the extraneous index_var columns, but that doesn't affect performance much).

I'm not calling pandas.tslib.array_to_timedelta64, to_dataframe is - the caller list is (sorry, I'm not sure of a better way to show this):

forcing_data.isel(lat=lat, lon=lon).values() returns a ValuesView, which scikit-learn doesn't like. However, forcing_data.isel(lat=lat, lon=lon).to_array().T seems to work..

@shoyer: Where does times come from in that code?

Example forcing data:

<xarray.Dataset> Dimensions: (lat: 360, lon: 720, time: 2928) Coordinates: * lon (lon) float64 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ... * lat (lat) float64 -89.75 -89.25 -88.75 -88.25 -87.75 -87.25 -86.75 ... * time (time) datetime64[ns] 2012-01-01 2012-01-01T03:00:00 ... Data variables: SWdown (time, lat, lon) float64 446.5 444.9 445.3 447.8 452.4 456.3 ...

Where there might be an arbitrary number of data variables, and the scikit-learn input would be time (rows) by data variables (columns). I'm currently doing this:

``` python def predict_gridded(model, forcing_data, flux_vars): """predict model results for gridded data

:model: TODO
:data: TODO
:returns: TODO

"""
# set prediction metadata
prediction = forcing_data[list(forcing_data.coords)]

# Arrays like (var, lon, lat, time)
result = np.full([len(flux_vars),
                  forcing_data.dims['lon'],
                  forcing_data.dims['lat'],
                  forcing_data.dims['time']],
                 np.nan)
print("predicting for lon: ")
for lon in range(len(forcing_data['lon'])):
    print(lon, end=', ')
    for lat in range(len(forcing_data['lat'])):
        result[:, lon, lat, :] = model.predict(
            forcing_data.isel(lat=lat, lon=lon)
                        .to_dataframe()
                        .drop(['lat', 'lon'], axis=1)
        ).T
print("")
for i, fv in enumerate(flux_vars):
    prediction.update(
        {fv: xr.DataArray(result[i, :, :, :], 
                          dims=['lon', 'lat', 'time'],
                          coords=forcing_data.coords)
        }
    )

return prediction

```

and I think it's working (still debugging, and it's pretty slow running)

I want to be able to run a scikit-learn model over a bunch of variables in a 3D (lat/lon/time) dataset, and return values for each coordinate point. Is something like this multi-dimensional groupby required (I'm thinking groupby(lat, lon) => 2D matrices that can be fed straight into scikit-learn), or is there already some other mechanism that could achieve something like this? Or is the best way at the moment just to create a null dataset, and loop over lat/lon and fill in the blanks as you go?

Not sure if it's worth a separate feature request, but it'd be great to have a drop option in .isel() as well. For example:

``` python

ds <xarray.Dataset> Dimensions: (lat: 360, lon: 720, time: 2928) Coordinates: * lon (lon) float64 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ... * lat (lat) float64 -89.75 -89.25 -88.75 -88.25 -87.75 -87.25 -86.75 ... * time (time) datetime64[ns] 2012-01-01 2012-01-01T03:00:00 ... Data variables: dswrf (time, lat, lon) float64 446.5 444.9 445.3 447.8 452.4 456.3 ...

ds.isel(lat=0, lon=0, drop=True) <xarray.Dataset> Dimensions: (time: 2928) Coordinates: * time (time) datetime64[ns] 2012-01-01 2012-01-01T03:00:00 ... Data variables: dswrf (time) float64 446.5 444.9 445.3 447.8 452.4 456.3 ... ```

This would be useful for making .to_dataframe() more efficient, too, I suspect, since I have to currently do forcing.isel(lat=0, lon=0).to_dataframe()[['dswrf']] to remove the (constant) lat/lon data from the dataframe.

Sorry, this is just something I stumbled across in a workshop, it's not a problem that's affecting my work, so I don't have much time for it. But I'll make the MOM developers aware of it - they might have an interest in helping out.

@shoyer: the file is publically available, have you tried checking it on your set-up? The file is about 35Mb..

Same problem with scipy - I'm not sure what I need to install to get pydap and h5netcdf working...

I updated to 1.2.1 from pypi, it still fails.

$ conda list |grep -i cdf libnetcdf 4.3.3.1 1 netcdf4 1.1.9 np110py34_0

xray.version.version == '0.6.1' from anaconda, on python 3.4

and?

The lat/x, long/y is a bit odd, yeah, I think the file is set up a bit wrong, but I'm trying to stick to the format I'm given.

Reference height, at least, isn't really a co-ordinate. It's the value of the tower measurements. The other data is really (implicitly) at elevation+screen height, I think.

But yeah, that's basically what I came up with, just copy the relevant coords for each variable from the old dataset.

Hrm, I guess I can just do

python new_ds = old_ds[GEO_VARS] for c in old_ds.coords: if c not in new_ds.coords: new_ds.coords[c] = old_ds.coords[c]

Still, it might be nice to have this built in.

Fair enough, I wasn't aware list(ds.dims) worked..

.keys() returns a keysview though, not a list:

``` python In [4]: ds.dims.keys() Out[4]: KeysView(Frozen(SortedKeysDict({'time': 70128, 'z': 1, 'x': 1, 'y': 1})))

In [5]: ds.dims.keys()[0]

TypeError Traceback (most recent call last) <ipython-input-5-874ed46e2639> in <module>() ----> 1 ds.dims.keys()[0]

TypeError: 'KeysView' object does not support indexing ```

Other envs aren't affected. I guess this is a conda problem, not an xray problem. Easy enough to re-create an env, I guess.

blegh... just noticed the different dates. Never mind :P

Thanks for the help again

Ah, cool, that looks good, however, I think there might be a bug somewhere. Here's what I was getting originally:

``` In [62]: new_data['time'].encoding

Out[62]: {}

In [60]: new_data.to_netcdf('data/tumba_site_mean_2_year.nc') ```

resulted in:

``` $ ncdump ../projects/synthetic_forcings/data/tumba_site_mean_2_year.nc|grep time --context=2 netcdf tumba_site_mean_2_year { dimensions: time = 35088 ; y = 1 ; x = 1 ; variables: ... float time(time) ; time:calendar = "proleptic_gregorian" ; time:units = "minutes since 2000-01-01 00:30:00" ; ...

time = 0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600, 630, 660, 690, 720, 750, 780, 810, 840, ```

And if I copy the encoding from the original file, as it was loaded:

``` In [65]: new_data['time'].encoding = data['time'].encoding

In [66]: new_data['time'].encoding

Out[66]: {'dtype': dtype('>f8'), 'units': 'seconds since 2002-01-01 00:30:00'}

In [67]: new_data.to_netcdf('data/tumba_site_mean_2_year.nc') ```

results in

``` $ ncdump ../projects/synthetic_forcings/data/tumba_site_mean_2_year.nc|grep time --context=1 dimensions: time = 35088 ; y = 1 ; -- variables: ... double time(time) ; time:calendar = "proleptic_gregorian" ; time:units = "seconds since 2002-01-01T00:30:00" ;

... time = -63158400, -63156600, -63154800, -63153000, -63151200, -63149400, -63147600, -63145800, -63144000, -63142200, -63140400, -63138600, ```

Now the units are right, but the values are way off. I can't see anything obvious missing from the encoding, compared to the xray docs, but I'm not sure how it works.

Also, since seconds are the base SI unit for time, I think it would be sensible to use seconds by default, if no encoding is given.

Ah, yep, making the dimension using data.coords['timeofday'] = ('time', [np.timedelta64(60 * int(h) + int(m), 'm') for h,m in zip(data['time.hour'], data['time.minute'])]) works. Thanks for all the help :)

Ok, weird. That example works for me, but even if I take a really short slice of my data set, the same thing won't work:

``` In [61]: d = data.sel(time=slice('2002-01-01','2002-01-03')) d

Out[61]:

<xray.Dataset> Dimensions: (time: 143, timeofday: 70128, x: 1, y: 1, z: 1) Coordinates: * x (x) >f8 1.0 * y (y) >f8 1.0 * z (z) >f8 1.0 * time (time) datetime64[ns] 2002-01-01T00:30:00 ... * timeofday (timeofday) timedelta64[ns] 1800000000000 nanoseconds ... Data variables: SWdown (time, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.58 ... Rainf_qc (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... SWdown_qc (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... Tair (time, z, y, x) float64 282.9 282.9 282.7 282.6 282.4 281.7 281.0 ... Tair_qc (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... LWdown (time, y, x) float64 296.7 297.3 297.3 297.3 297.2 295.9 294.5 ... PSurf_qc (time, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... latitude (y, x) float64 -35.66 Wind (time, z, y, x) float64 2.2 2.188 1.9 2.2 2.5 2.5 2.5 2.25 2.0 2.35 ... LWdown_qc (time, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... Rainf (time, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... Qair_qc (time, y, x) float64 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... longitude (y, x) float64 148.2 PSurf (time, y, x) float64 8.783e+04 8.783e+04 8.782e+04 8.781e+04 ... reference_height (y, x) float64 70.0 elevation (y, x) float64 1.2e+03 Qair (time, z, y, x) float64 0.00448 0.004608 0.004692 0.004781 ... Wind_qc (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... Attributes: Production_time: 2012-09-27 12:44:42 Production_source: PALS automated netcdf conversion Contact: palshelp@gmail.com PALS_fluxtower_template_version: 1.0.2 PALS_dataset_name: TumbaFluxnet PALS_dataset_version: 1.4

In [62]: d.groupby('timeofday').mean('time') ```

That last command will not complete - it will run for minutes. Not really sure how to debug that behaviour.

Perhaps it's to do with the long/lat/height variables that really should be coordinates (I'm just using the data as it came, but I can clean that, if necessary)

same problem with numpy.timedelta64 too.

Dunno if this is related to the ds['time.time'] problem, but I tried creating the daily_cycle using a pandas.Timedelta as the index (timeofday), but it also appeared to just hang indefinitely when doing the data.groupby('timeofday').mean('time') call..

Nice.

Ok, I have hit a stumbling block, and this is much more of a support request, so feel free to direct me else where, but since we're on the topic, I want to do something like:

start = 2002 n_years = 4 new_data = [] for year in range(start, start + n_years): days = 365 if year%4 else 365 for d in range(days): day_data = mean + annual_cycle.isel(dayofyear=d) + daily_cycle day_data.coords['time'] = datetime.datetime(year,1,1) + datetime.timedelta(day=d, hour=day_data.timeofday.hour, minute=day_data.timeofday.minute) new_data.append(day_data) xray.concat(new_data)

where mean, annual_cycle, and daily_cycle are overall mean, annual cycle at daily resolution, and daily cycle at 30 minute resolution (the latter two bias corrected by subtracting the mean). I'm trying to make a synthetic dataset 4 years long that only includes the mean, seasonal, and daily cycles, but no other variability.

The assignment of day_data['time'] fails because the day_data.timeofday.hour (and .minute) don't work. These are datetime.times, is there an efficient way of converting these to datetime.timedeltas, without first manually taking them out of the DataArray?

Ah, cool, thanks for that link, I missed that in the docs.

One thing that would be nice (in both pandas and xray) is a time.timeofday. I can't figure out how to do it with time.hour and time.minute - I need half-hourly resolution averaging. time.time does something in xray, but it seems to never complete, and it doesn't work at all in pandas.

Heh, I meant the pandas docs - they don't specify the rule argument format either

time.month and time.hour do exactly what I need. They aren't mentioned in the docs at http://xray.readthedocs.org/en/stable/groupby.html, and I'm not sure how I'd guess that they exist, so perhaps they should be added to that page? It doesn't appear to be something that exists in pandas..

Unfortunately I'm not familiar enough with pd.resample and pd.TeimGrouper to know the difference in what they can do. resample looks like it would cover my use-cases, although the docs are pretty limited, and don't actually specify the format of the rule format...

One thing that I would like to be able to do that is not covered by resample, and might be covered by TimeGrouper is to group over month only (not month and year), in order to create a plot of mean seasonal cycle (at monthly resolution), or similarly, a daily cycle at hourly resolution. I haven't figured out if I can do that with TimeGrouper yet though.

Looks good to me. I don't know enough to be able to comment on the API question.