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