In [ ]:
import pyproj
crs = pyproj.CRS.from_user_input('epsg:3003')
crs
Out[ ]:
<Projected CRS: EPSG:3003>
Name: Monte Mario / Italy zone 1
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: Italy - onshore and offshore - west of 12°E.
- bounds: (5.93, 36.53, 12.0, 47.04)
Coordinate Operation:
- name: Italy zone 1
- method: Transverse Mercator
Datum: Monte Mario
- Ellipsoid: International 1924
- Prime Meridian: Greenwich
In [ ]:
# Create a pyproj CRS object
crs = pyproj.CRS.from_user_input('epsg:2154') # lambert conformal
crs = pyproj.CRS.from_user_input('epsg:3003') # monte mario

# Extract projection parameters from the pyproj CRS object
projection_name = crs.name
geographic_crs_name = crs.geodetic_crs.name
horizontal_datum_name = crs.datum.name
reference_ellipsoid_name = crs.ellipsoid.name
prime_meridian_name = crs.prime_meridian.name
false_easting = crs.coordinate_operation.params[0].value
false_northing = crs.coordinate_operation.params[1].value
central_meridian = crs.coordinate_operation.params[2].value
scale_factor = crs.coordinate_operation.params[3].value
latitude_of_origin = crs.coordinate_operation.params[4].value
semi_major_axis = crs.ellipsoid.semi_major_metre
inverse_flattening = crs.ellipsoid.inverse_flattening

# Create a dictionary to hold CF Grid Mapping Attributes
cf_attributes = {
    'grid_mapping_name': 'transverse_mercator',
    'projected_crs_name': projection_name,
    'geographic_crs_name': geographic_crs_name,
    'horizontal_datum_name': horizontal_datum_name,
    'reference_ellipsoid_name': reference_ellipsoid_name,
    'prime_meridian_name': prime_meridian_name,
    'false_easting': false_easting,
    'false_northing': false_northing,
    'longitude_of_central_meridian': central_meridian,
    'scale_factor_at_central_meridian': scale_factor,
    'latitude_of_projection_origin': latitude_of_origin,
    'semi_major_axis': semi_major_axis,
    'inverse_flattening': inverse_flattening,
    # Add other CF attributes here as needed based on the CSV files
}

# Print the CF Grid Mapping Attributes
cf_attributes
Out[ ]:
{'grid_mapping_name': 'transverse_mercator',
 'projected_crs_name': 'Monte Mario / Italy zone 1',
 'geographic_crs_name': 'Monte Mario',
 'horizontal_datum_name': 'Monte Mario',
 'reference_ellipsoid_name': 'International 1924',
 'prime_meridian_name': 'Greenwich',
 'false_easting': 0.0,
 'false_northing': 9.0,
 'longitude_of_central_meridian': 0.9996,
 'scale_factor_at_central_meridian': 1500000.0,
 'latitude_of_projection_origin': 0.0,
 'semi_major_axis': 6378388.0,
 'inverse_flattening': 297.0}
In [ ]:
import numpy as np
x, y = np.meshgrid(np.linspace(4770000 ,4770000+100000,100), np.linspace(1620000,1620000+100000,100))
In [ ]:
import xarray as xr
import numpy as np

# Create a xarray Dataset
# Define variables

ds = xr.Dataset({
    'x': xr.DataArray(x, dims=['x','y'], attrs={
        'units': "m",
        'long_name': "x coordinate of projection",
        'standard_name': "projection_x_coordinate"
    }), 
    'y': xr.DataArray(y, dims=['x','y'], attrs={
        'units': "m",
        'long_name': "y coordinate of projection",
        'standard_name': "projection_y_coordinate"
    })
})
ds.attrs = cf_attributes

# Save the xarray Dataset to a NetCDF file
ds
Out[ ]:
<xarray.Dataset> Size: 160kB
Dimensions:  (x: 100, y: 100)
Coordinates:
    x        (x, y) float64 80kB 4.77e+06 4.771e+06 ... 4.869e+06 4.87e+06
    y        (x, y) float64 80kB 1.62e+06 1.62e+06 ... 1.72e+06 1.72e+06
Data variables:
    *empty*
Attributes: (12/13)
    grid_mapping_name:                 transverse_mercator
    projected_crs_name:                Monte Mario / Italy zone 1
    geographic_crs_name:               Monte Mario
    horizontal_datum_name:             Monte Mario
    reference_ellipsoid_name:          International 1924
    prime_meridian_name:               Greenwich
    ...                                ...
    false_northing:                    9.0
    longitude_of_central_meridian:     0.9996
    scale_factor_at_central_meridian:  1500000.0
    latitude_of_projection_origin:     0.0
    semi_major_axis:                   6378388.0
    inverse_flattening:                297.0
xarray.Dataset
    • x: 100
    • y: 100
    • x
      (x, y)
      float64
      4.77e+06 4.771e+06 ... 4.87e+06
      units :
      m
      long_name :
      x coordinate of projection
      standard_name :
      projection_x_coordinate
      array([[4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ],
             [4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ],
             [4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ],
             ...,
             [4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ],
             [4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ],
             [4770000.       , 4771010.1010101, 4772020.2020202, ...,
              4867979.7979798, 4868989.8989899, 4870000.       ]])
    • y
      (x, y)
      float64
      1.62e+06 1.62e+06 ... 1.72e+06
      units :
      m
      long_name :
      y coordinate of projection
      standard_name :
      projection_y_coordinate
      array([[1620000.       , 1620000.       , 1620000.       , ...,
              1620000.       , 1620000.       , 1620000.       ],
             [1621010.1010101, 1621010.1010101, 1621010.1010101, ...,
              1621010.1010101, 1621010.1010101, 1621010.1010101],
             [1622020.2020202, 1622020.2020202, 1622020.2020202, ...,
              1622020.2020202, 1622020.2020202, 1622020.2020202],
             ...,
             [1717979.7979798, 1717979.7979798, 1717979.7979798, ...,
              1717979.7979798, 1717979.7979798, 1717979.7979798],
             [1718989.8989899, 1718989.8989899, 1718989.8989899, ...,
              1718989.8989899, 1718989.8989899, 1718989.8989899],
             [1720000.       , 1720000.       , 1720000.       , ...,
              1720000.       , 1720000.       , 1720000.       ]])
      • grid_mapping_name :
        transverse_mercator
        projected_crs_name :
        Monte Mario / Italy zone 1
        geographic_crs_name :
        Monte Mario
        horizontal_datum_name :
        Monte Mario
        reference_ellipsoid_name :
        International 1924
        prime_meridian_name :
        Greenwich
        false_easting :
        0.0
        false_northing :
        9.0
        longitude_of_central_meridian :
        0.9996
        scale_factor_at_central_meridian :
        1500000.0
        latitude_of_projection_origin :
        0.0
        semi_major_axis :
        6378388.0
        inverse_flattening :
        297.0
      In [ ]:
      transformer = pyproj.Transformer.from_crs(crs, 'epsg:4326')
      
      In [ ]:
      import geoviews as gv
      import geopandas as gp
      import hvplot.pandas
      gv.extension('bokeh')
      lon, lat  = transformer.transform(ds.x, ds.y)
      
      map = gp.read_file(gp.datasets.get_path('naturalearth_lowres'))
      map_ = map.hvplot(geo=True)
      (map_ * gv.Points((lon.ravel(), lat.ravel()) ).opts(color = 'r')).opts(width=600, height=600)
      
      /tmp/ipykernel_26085/1761570461.py:7: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
        map = gp.read_file(gp.datasets.get_path('naturalearth_lowres'))
      
      Out[ ]: