Changing EPSG values from .shp file

Specific Domains Geo
1

I am trying to convert my .shp file from an EPSG value of 26917 to a EPSG value of 4326. These seems to be an issue with my .shp file not having the EPSG built in. I have tried changing the EPSG using the code below but I receive this error: No CRS found in the GeoDataFrame. Setting to EPSG:4326.
Cannot reproject geometries because no source CRS is defined.
And when I try to plot the geometries later the EPSG value never changes. Any suggestions or advice would be greatly appreciated.

using ArchGDAL
using GeoDataFrames

Load the shapefile into a GeoDataFrame

shapefile_path = “…\Data & Resources for Ivy\Williford_10row_30m_rectangle.shp\Williford_10row_30m_rectangle.shp”
geo_df = GeoDataFrames.read(shapefile_path)

Check if the CRS exists in the GeoDataFrame

if in(:crs, names(geo_df))
println("Original CRS: ", geo_df.crs)
else
println(“No CRS found in the GeoDataFrame. Setting to EPSG:4326.”)
# Create a new CRS reference
crs = ArchGDAL.importEPSG(4326) # Adjust EPSG code as necessary
# Assign the CRS to the GeoDataFrame
geo_df[!, :crs] .= [crs] # Use broadcasting to assign CRS as a vector
end

Define the target CRS (e.g., WGS 84)

target_crs = ArchGDAL.importEPSG(4326) # Change to the desired target CRS if different

Extract geometries from the GeoDataFrame

geometries = geo_df.geom

Reproject the geometries if the source CRS is set

if in(:crs, names(geo_df))
source_crs = geo_df.crs
# Reproject the geometries
reprojected_geometries = ArchGDAL.reproject(geometries, source_crs, target_crs)

# Assign the reprojected geometries back to the GeoDataFrame
geo_df.geom .= reprojected_geometries  # Use broadcasting to assign

# Output the WKT of the reprojected geometries to verify
println("Reprojected Geometries:")
for geom in geo_df.geom
    wkt_geom = ArchGDAL.toWKT(geom)
    println(wkt_geom)
end

else
println(“Cannot reproject geometries because no source CRS is defined.”)
end

2

The .shp file doesn’t store the CRS. You need to find the .proj file in the same directory. If you don’t have one, then you need to ask the author of the file for the CRS.

Also, take a look at

https://juliaearth.github.io/geospatial-data-science-with-julia/06-projections.html

3

If you know the source EPSG just specify it manually as e.g. EPSG(1234). You may need to do using GeoFormatTypes for that to work.

4

What about something like this [NOT TESTED]?

using GeoDataFrames
import GeoDataFrames as GFT
using Proj

path_2_shapefile = "/example/path/my_shapefile.shp"

# assign EPSG code to the shapefile
gdf = GeoDataFrames.read(path_2_shapefile)
GeoDataFrames.write(path_2_shapefile, gdf; crs=GFT.EPSG(26917))

source_crs1 = GFT.EPSG(26917)
target_crs1 = GFT.EPSG(4326)
trans = Proj.Transformation(source_crs1, target_crs1; always_xy=true)

gdf[:, :geometry] = GO.transform.(trans, gdf.geometry)

path_2_shapefile = "/example/path/my_shapefile_4326.shp"
GeoDataFrames.write(path_2_shapefile, gdf; crs=GFT.EPSG(4326))
1 Like
5

GO.reproject might simplify that too, as you don’t need to define the transformation.

(I was trying to stick to ops use of ArchGDAL.jl, but as author of both reproject methods and transform I would use GO.reproject now)

Also needs import GeoFormatTypes as GFT, small typo there in the package name

6

Hi @ivyc, welcome to the community. May I ask where you got the idea to store the crs as a column in the DataFrame?

The CRS in a (Geo)DataFrame can be stored in its metadata and/or in the stored geometries. It can also be hardcoded during writing (crs keyword). I just made a new release of GeoDataFrames (v0.3.11), which enables the following:

using GeoDataFrames  # no other imports needed

gdf = GeoDataFrames.read("test.shp")
crs(gdf)  # nothing, because no .prj was found
df = reproject(gdf, EPSG(26917), EPSG(4326))
crs(df)  # EPSG{1}((4326,))

GeoDataFrames.write("test_reprojected.shp", df)

As pointed out by Júlio, shapefiles requires a separate file for storing their crs. I would suggest to move away from shapefiles, and use something like a geopackage. Just GeoDataFrames.write("test.gpkg", df) would work for that.

Let me know if this works for you, and whether you have any other questions.

1 Like
7

To add another option, this how we do that with GMT.jl The test file is a tinny one from the GDAL autotests, renamed to .shp.zip to simplify things. But since .zip files (or .shp) are not accepted here I cannot provide that data file

D = gmtread("poly.shp.zip")
Vector{GMTdataset} with 10 segments
...
PROJ: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
...

Now reproject it to geogs

xy2lonlat(D, t_srs=4326)
Vector{GMTdataset} with 10 segments
...
PROJ: +proj=longlat +datum=WGS84 +no_defs

20×2 GMTdataset{Float64, 2}
 Row │     lon      lat
─────┼──────────────────
   1 │ 162.859  87.2616
   2 │ 162.887  87.2613
   3 │ 162.9    87.2604
...