andypbarrett
commented
1 year ago Hi
1. Why the latitude values aren't negative (is the function just for Northern Hem? )
The function converts projected coordinates (x, y) for both Northern and Southern Polar Stereographic projections. I tested the Southern Hemisphere case for the projected grid origin (The South Pole) for the S. Hem. (x=0., y=0.), which is (lon=0., lat=-90.)
>>> from polar_convert.constants import NORTH, SOUTH
>>> from polar_convert import polar_xy_to_lonlat
>>> x, y = 0, 0
>>> true_scale_lat = 70
>>> re = 6378.137
>>> e = 0.01671
>>> hemisphere = NORTH
>>> polar_xy_to_lonlat(x, y, true_scale_lat, re, e, hemisphere)
[0.0, 90.0]
>>> hemisphere = SOUTH
>>> true_scale_lat = -70 # For the S. Hem, the latitude of true scale is -70 N
>>> polar_xy_to_lonlat(x, y, true_scale_lat, re, e, hemisphere)
[0.0, -90.0]However, polar_xy_to_lonlat expects projected coordinates in kilometers. The x,y coordinates in the dataset you are using appear to be in meters. I cannot see the units attributes for the x and y coordinates, so you should check this. Even with hemisphere=SOUTH, the first elements of the x and y coordinates, without conversion to km, gives a location close to the North Pole.
>>> x = -5.219e+05
>>> y = -1.891e+06
>>> polar_xy_to_lonlat(x, y, true_scale_lat, re, e, SOUTH)
[195.42905457549008, 89.27762998161299]One way to check if this is what is happening is to find the minimum latitude in the lat grid you have calculated.
ds.lat.min()Another check is the plot the latitude grid, either over a map or just as an image to check both the range of latitude values, and the spatial extent of the points.
Once x and y are converted to km, you get the result you expect for a S. Hem grid.
>>> y = -1.891e+06 * 1e-3
>>> x = -5.219e+05 * 1e-3
>>> polar_xy_to_lonlat(x, y, true_scale_lat, re, e, hemisphere)
[195.42905457549008, -71.98213283556814]However, I do note that -71.9 degrees is a much higher latitude that our standard polar stereographic grids. See here: https://nsidc.org/data/user-resources/help-center/guide-nsidcs-polar-stereographic-projection
I suggest confirming that the grid is a polar stereographic grid. Many of our AMSR2 products are EASE2 (https://nsidc.org/data/user-resources/help-center/guide-ease-grids).
2. Why am I not able to modify the data using these lat-lon?
The xarray.DataArray.sel() method requires that indexers are vectors. lat is an array (2D). Even though you have assigned lat and lon as coordinate variables, only x, y and time can be used for indexing because they are the only coordinates that are (1D) vectors.
Which approach you use to subset the data, depends on what you want to do with the subsetted data.
If you want to extract a rectangular region of data, you can determine the lat, lon coordinates of the corner points of that region, and then use polar_lonlat_to_xy to get a bounding box in projected coordinates . When I do this, I usually use a combination of Google Maps or Google Earth, plotting the corner points over the data using matplotlib and cartopy, and trial and error. You could also use something like shapely or geopandas to define a polygon, convert that polygon to xy coordinates, and then use the bounds methods to get the best-fit bounding-box for that polygon.
This will yield a bounding-box, similar to [minx, miny, maxx, maxy]. These bounds can be used to subset the data using .sel(). For example:
ds_subset = ds.Z.sel(x=slice(minx, maxx), y=slice(miny, maxy))However, if you want to get a latitude band to calculate a zonal mean, for example, you can mask the latitude band you want and calculate the mean.
ds_masked = ds.Z.where((ds.lat >= 82) & (ls.lat < 83))
zonal_mean = ds_masked.mean()Alternatively, you can reproject and resample the polar stereographic grid to a latitude-longitude grid. polar_ij_to_lonlat can be used to reproject the grid points to latitude-longitude. These would then need to be interpolated to a regular grid. You could also use rioxarray to do this.
Note, that if you use rioxarray or some tool other than the ones here, and if the data are in either Polar Stereographic or EASE-Grid projection that use the Hugh 1980 or International 1924 Authalic Sphere as their datums, you will need to specify the CRS as a proj4 string or WKT, not as an EPSG code. This is because these two datums have been depreciated by EPSG and reprojection software such as gdal silently use projection definitions that have WGS84 as the datum. This yields incorrect results. We are working with EPSG to solve this issue.
gkb999
Hi there, I have been working with NSIDC datasets for ages now, but came across this a while ago. Good one firstly.
Second things, does the cocnversion work for AMSR2 data too? As they too follwo the similar NSIDC structure.
My aim: Add additional dimensions as lat-lon along with x and y for the dataset (nc file). Waht I did:
What I got:
So, looks all fine, but 1) the latitudes aren't negative why? 2) i'm trying to assign the obtained arrays to my original dataset.
What I tried:
Looks fine:
My dataset now looks like:
But, when I tried
I get the follwoing error:
So I'm expecting 2 things to be sorted. 1) Why the latitude values aren't negative (is the function just for Northern Hem? ) 2) Why am I not able to modify the data using these lat-lon?