Coordinate priority with non-spherical CRS

[bekozi]

@mkmitchell reported an issue with a NARR dataset having a Lambert Conformal Conic projection. The lat/lon coordinate variables were marked with axis x/y attributes. The Cartesian x/y coordinate variables had no axis attributes. ocgis errored out during a subset as it chose the LCC projection correctly but picked the wrong coordinates because of axis attribute priority during coordinate selection. This was fixed by updating the dimension map prior to the subset.

To fix this, coordinate variable selection should attempt to match coordinate variable units against the units recommended by the coordinate system prior to using the axis attributes.

[mkmitchell]

I was able to bypass the error by using your suggestion

rdc = RequestDataset(uri, 'variable')
rdc.dimension_map.set_variable('x', 'x')
rdc.dimension_map.set_variable('y', 'y')

but it only worked for a bit. I'm getting the same error again.

ocgis.exc.ExtentError: A subset operation on variable "None" returned empty. This typically means the selection geometry falls outside the spatial domain of the target dataset.

I'm attempting to sum the variable by day/month/year for every shape in the shapefile using the following command:

ops = OcgOperations(dataset=rdc, output_format='shp', time_region={'month': [1,2,3,4,9,10,11,12]}, spatial_operation='clip', geom=shp, calc=calc,
                calc_raw=True, aggregate=True, calc_grouping=['day', 'month', 'year'], prefix=prefix)

Thanks in advance!

[bekozi]

Bummer! A few quick questions:

• Did the netcdf metadata change?
• Are you using a different shapefile?
• Did you update ocgis recently?

I cannot think of any recent changes that would cause this. It's possible (if the shapfile changed) that you do have some empty subsets in there. In that case, you can try setting allow_empty=True in the operations call.

[mkmitchell]

• I've created multiple NetCDF files but have also used original NARR data with the same result.

It must be the shapefile. I get the error I posted about with one. It succeeds with a few others.
I've also seen ValueError: Record's geometry type does not match collection schema's geometry type: 'Point' != 'MultiPoint' with a different one.

[mkmitchell]

It was definitely the shapefile. It appears very small features were giving me the error.

[bekozi]

Thanks for the update. Sorry for the slow reply - was taking some PTO.

Do you think this is a problem that needs to be fixed in ocgis? I could take a look at the shapefile if that would help.

[mkmitchell]

Take your much deserved PTO!

That's completely up to you. I ended up iterating through selections of my shapefile. If it worked I selected more and if it didn't I deselected a few I previously added.

It may be nice to have ocgis keep running the calculation and report a list of features that didn't process because of whatever exception.

[bekozi]

Done! :wink: Your approach makes sense. When you get a chance, could you pass along the shp and netcdf (or just the spatial resolution of the grid if that's easier)? What you described would be a nice feature.