import geopandas as gpd

# Filepath to the addresses Shapefile
fp = r"/home/geo/addresses.shp"
fp = r"D:\KOODIT\Opetus\Automating-GIS-processes\AutoGIS-Sphinx\source\data\addresses.shp"

# Read data
data = gpd.read_file(fp)

"""
Explanation: Projections: Converting from projection to another
A map projection is a systematic transformation of the latitudes and longitudes into a plain surface. As map projections of gis-layers are fairly often defined differently (i.e. they do not match), it is a common procedure to redefine the map projections to be identical in both layers. It is important that the layers have the same projection as it makes it possible to analyze the spatial relationships between layer, such as conduct the Point in Polygon spatial query (which we will try next). 
Changing the coordinate reference system of a layer in Geopandas
Defining a projection and changing it is easy in Geopandas. Let's continue working with our  address points, and change the Coordinate Reference System (CRS) from WGS84 into a projection called ETRS GK-25 (EPSG:3879) which uses a Gauss-Krüger projection that is (sometimes) used in Finland. 

Let's first read the data from the Shapefile that we created previously
End of explanation
"""


data.crs

"""
Explanation: Let's check what is the current CRS of our layer
End of explanation
"""


data['geometry'].head()

"""
Explanation: Okey, so it is WGS84 (i.e. EPSG: 4326). 

Let's also check the values in our geometry column
End of explanation
"""


# Let's take a copy of our layer
data_proj = data.copy()

# Reproject the geometries by replacing the values with projected ones
data_proj['geometry'] = data_proj['geometry'].to_crs(epsg=3879)

"""
Explanation: Okey, so they indeed look like lat-lon values. 


Let's convert those geometries into ETRS GK-25 projection (EPSG: 3879). Changing the projection is really easy to do in Geopandas with .to_crs() -function. As an input for the function, you should define the column containing the geometries, i.e. geometry in this case, and a epgs value of the projection that you want to use. 


Note: there is also possibility to pass the projection information as proj4 strings or dictionaries, see more here
End of explanation
"""


data_proj['geometry'].head()

"""
Explanation: Let's see how they look now
End of explanation
"""


# Plot the WGS84 
data.plot()

# Plot the one with ETRS GK-25 projection
data_proj.plot()

"""
Explanation: And here we go, the numbers have changed! Now we have successfully changed the projection of our layer into a new one. 

Let's still compare the layers visually
End of explanation
"""


from fiona.crs import from_epsg

# Determine the CRS of the GeoDataFrame
data_proj.crs = from_epsg(3879)

# Let's see what we have
print(data_proj.crs)

"""
Explanation: Indeed, they look different and our re-projected one looks much better in Finland (not so stretced as in WGS84). 

Now we still need to change the crs of our GeoDataFrame into EPSG 3879 as now we only modified the values of the geometry column. We can take use of fiona's from_epsg -function.
End of explanation
"""


# Pass the coordinate information
data_proj.crs = {'y_0': 0, 'no_defs': True, 'x_0': 25500000, 'k': 1, 'lat_0': 0, 'units': 'm', 'lon_0': 25, 'ellps': 'GRS80', 'proj': 'tmerc'}

# Check it
print(data_proj.crs)

"""
Explanation: Notice: The above works for most EPSG codes but ETRS GK-25 projection is rather rare, so we still need to specify and make sure that .prj file is having correct coordinate system information by passing a proj4 dictionary below into it (otherwise the .prj file would be empty):
End of explanation
"""


# Ouput file path
outfp = r"/home/geo/addresses_epsg3879.shp"

# Save to disk
data_proj.to_file(outfp)

"""
Explanation: Finally, let's save our projected layer into a Shapefile so that we can use it later.
End of explanation
"""


