I want to find the corresponding ice concentration values that are present at a series of lat and lon coordinates (such as lon=-74.7732772827148 and lat=-69.3250350952148). However, the dataset I am using (https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-ice-concentration?tab=form) has some of the following attributes:
double xc(xc) ;
xc:units = "km" ;
xc:long_name = "x coordinate of projection (eastings)" ;
xc:standard_name = "projection_x_coordinate" ;
double yc(yc) ;
yc:units = "km" ;
yc:long_name = "y coordinate of projection (northings)" ;
yc:standard_name = "projection_y_coordinate" ;
float lat(yc, xc) ;
lat:units = "degrees_north" ;
lat:long_name = "latitude coordinate" ;
lat:standard_name = "latitude" ;
float lon(yc, xc) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude coordinate" ;
lon:standard_name = "longitude" ;
int ice_conc(time, yc, xc) ;
ice_conc:_FillValue = -32767 ;
ice_conc:long_name = "fully filtered concentration of sea ice using atmospheric correction of brightness temperatures and open water filters" ;
ice_conc:standard_name = "sea_ice_area_fraction" ;
ice_conc:units = "%" ;
ice_conc:valid_min = 0 ;
ice_conc:valid_max = 10000 ;
ice_conc:grid_mapping = "Lambert_Azimuthal_Grid" ;
ice_conc:coordinates = "time lat lon" ;
ice_conc:ancillary_variables = "total_standard_error status_flag" ;
ice_conc:scale_factor = 0.01 ;
ice_conc:comment = "this field is the primary sea ice concentration estimate for this climate data record"
How do I find the corresponding x and y if lat and lon are both a function of x and y?
I have tried changing the projection using pyproj but it doesn't project correctly. As an example here is the following:
lambert_aea1 = {'proj': 'laea',
'lat_0':-90,
'lon_0':0,
'ellps': 'WGS84',
'datum': 'WGS84'}
inProj = Proj(init = 'epsg:4326')
outProj1 = Proj(lambert_aea1)
x1,y1 = np.array(transform(inProj,outProj1,lat_63_1[j],lon_63_1[j]))
the results are x1 =-1586154.039780751 and y1=598565.7767689897
but x and y values are in the range of -5387.5 and 5387.5
Is there another way to directly get the values rather than changing the projection?
If you look at the metadata for x and y in your file, the units are kilometers. For most projections, the units are either meters or feet.
If you multiply your x and y by 1,000 and compare the results, it is within the range you would expect.