Facing difficultyin converting Coordinate system from x,y to lat,lon

Retrievals
1

I am trying to convert x and y coordinates into latitude and longitude. Can someone here help me in guiding?

@kmuehlbauer @mgrover1

Below is the description of the metadata.

Thanks in advance.
Sandeep

dataset_2017
Out[39]: 
<xarray.Dataset>
Dimensions:      (time: 8760, y: 254, x: 256)
Coordinates:
  * time         (time) datetime64[ns] 2017-01-01T00:50:00 ... 2017-12-31T23:...
  * y            (y) float64 -4.171e+06 -4.17e+06 ... -3.919e+06 -3.918e+06
  * x            (x) float64 8.904e+04 9.004e+04 9.104e+04 ... 3.43e+05 3.44e+05
Data variables:
    spatial_ref  int32 ...
    RW           (time, y, x) float32 ...
Attributes:
    radarid:            10000
    formatversion:      3
    radolanversion:     2.18.3
    radarlocations:     ['boo', 'ros', 'emd', 'hnr', 'umd', 'pro', 'ess', 'fl...
    moduleflag:         1
    reanalysisversion:  2017.002
    intervalunit:       0

dataset_2017['spatial_ref'].attrs
Out[43]: 
{'crs_wkt': 'PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.314245179,
 'inverse_flattening': 298.257223563,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'unknown',
 'horizontal_datum_name': 'Unknown based on WGS 84 ellipsoid',
 'projected_crs_name': 'unknown',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 0.0,
 'false_northing': 0.0,
 'spatial_ref': 'PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]',
 'GeoTransform': '88537.83307814406 1000.0 0.0 -4171644.7242655726 0.0 1000.0'}


dataset_2017['spatial_ref'].attrs['crs_wkt']
Out[44]: 'PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]'
2

Hey

I am not sure this works in your case, but you could try direct georeferencing with the xarray Dataset:

proj = wradlib.georef.epsg_to_osr(4326) # Define projection | EPSG: 4326 for WGS84
dataset_2017_georef = dataset_2017.copy().pipe(wradlib.georef.georeference_dataset, proj=proj) # Georeference the data internally (link)

Or if this doesn’t work then georeference the dataset explicitly using the x and y data from your Dataset:

proj_source = wradlib.georef.projection.wkt_to_osr(dataset_2017[‘spatial_ref’].attrs[‘crs_wkt’]) # link
coords = wradlib.georef.projection.reproject(dataset_2017.x.values, dataset_2017.y.values, projection_source=proj_source, proj_target=proj) # From the polar stereographic to WGS84 (link, also example here)

Something like that. Then use the newly reprojected coordinates to plot or any other further steps you need to do.

Hope this helps!

Best regards
Jorma Rahu

1 Like
3

Dear Jorma Rahu,

Thank you for the reply.

I tried below code:

proj = wrl.georef.epsg_to_osr(4326)
proj_source = wrl.georef.projection.wkt_to_osr(dataset_2017['spatial_ref'].attrs['crs_wkt'])
coords = wrl.georef.projection.reproject(dataset_2017.x.values, dataset_2017.y.values, projection_source=proj_source, proj_target=proj)


coords = wrl.georef.projection.reproject(dataset_2017.x.values, dataset_2017.y.values, projection_source=proj_source, proj_target=proj)
Traceback (most recent call last):

  Cell In[15], line 1
    coords = wrl.georef.projection.reproject(dataset_2017.x.values, dataset_2017.y.values, projection_source=proj_source, proj_target=proj)

  File ~\anaconda3\envs\wradlib\Lib\site-packages\wradlib\georef\projection.py:288 in reproject
    raise TypeError("Incompatible X, Y inputs to 'reproject'")

TypeError: Incompatible X, Y inputs to 'reproject'
4

@SandeepAllampalli The error you are seeing is triggered when the shapes of X and Y arrays mismatch. Please check the shapes of the input arrays.

We could do better error messages.

5

A bit more context. You have X,Y as 1Dim. But when you want to project to latlon this will need 2D coordinates. So you would need to create a meshgrid from your 1d data and stack it (X,Y,2).

I can create an example along your data the next days, if my sparse explanations (sorry, late in the evening here) do not really help here.

1 Like
6

Hi @kmuehlbauer ,

dataset looks like this:

dataset_2017
Out[16]: 
<xarray.Dataset>
Dimensions:      (time: 8760, y: 254, x: 256)
Coordinates:
  * time         (time) datetime64[ns] 2017-01-01T00:50:00 ... 2017-12-31T23:...
  * y            (y) float64 -4.171e+06 -4.17e+06 ... -3.919e+06 -3.918e+06
  * x            (x) float64 8.904e+04 9.004e+04 9.104e+04 ... 3.43e+05 3.44e+05
Data variables:
    spatial_ref  int32 ...
    RW           (time, y, x) float32 ...
Attributes:
    radarid:            10000
    formatversion:      3
    radolanversion:     2.18.3
    radarlocations:     ['boo', 'ros', 'emd', 'hnr', 'umd', 'pro', 'ess', 'fl...
    moduleflag:         1
    reanalysisversion:  2017.002
    intervalunit:       0

dataset_2017['spatial_ref'].attrs
Out[43]: 
{'crs_wkt': 'PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.314245179,
 'inverse_flattening': 298.257223563,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'unknown',
 'horizontal_datum_name': 'Unknown based on WGS 84 ellipsoid',
 'projected_crs_name': 'unknown',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 0.0,
 'false_northing': 0.0,
 'spatial_ref': 'PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",10],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]',
 'GeoTransform': '88537.83307814406 1000.0 0.0 -4171644.7242655726 0.0 1000.0'}
7

So according to Kai’s answer you just need to meshgrid the coords:

proj = wrl.georef.epsg_to_osr(4326)
proj_source = wrl.georef.projection.wkt_to_osr(dataset_2017[‘spatial_ref’].attrs[‘crs_wkt’])

xx, yy = np.meshgrid(dataset_2017.x.values, dataset_2017.y.values)
coords = wrl.georef.projection.reproject(xx, yy, projection_source=proj_source, proj_target=proj)

This did not produce error when I tested it.

2 Likes
8

Thank you @jorahu , it worked.

2 Likes