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.0In [ ]:
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[Â ]: