[FIXED] Finding the xy values from lat lon from SSMI/S sea ice concentration data

Issue

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?

Solution

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.

Answered By - snowman2