Get adequate projected coordinate system from geographic coordinate system across globe with GeoPandas in Python

coordinate systemgeopandaspython

The Situation

I have a global multiband-raster of 10×10 km resolution at the equator. I need to turn this raster into a GeoDataframe for further processing (where every band would become one column of values). This GeoDataframe is a global grid, also 10×10 km resolution at the equator. I calculate the centroid of every cell/geometry of the GeoDataframe and retrieve the underlying value of the raster with the point_query function from shapely.

The Problem

The global raster is in the geographic coordinate system EPSG:4326. Employing the shapely centroid function is imprecise when used with a geographic coordinate system, mostly close to the poles (if I understood correctly). For the imprecision, I am referring to the following UserWarning raised by GeoPandas:

UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

The Question

Is there a way to assign the correct projected coordinate reference system (e.g., UTM) globally to every "group" of grid cells, and therefore not get the GeoPandas UserWarning raised?

What I have tried

The code below retrieves the centroid of every country in latitude and longitude (could also be any arbitrary grid, e.g. of 10×10 degrees), and generates a UTM code. The problem is, that the calculation of this centroid is already imprecise.

        # Take the UTM coordinate system to correctly calculate the area, since epsg4326 is not a planar crs. See https://gis.stackexchange.com/questions/365584/convert-utm-zone-into-epsg-code
        # set coordinate system for area calculation for points
        country_center_x = TM_border.get_borders().centroid.iloc[
            0].x  # get the x coordinate of the center of the country
        country_center_y = TM_border.get_borders().centroid.iloc[
            0].y  # get the y coordinate of the center of country
        specific_utm = utm.from_latlon(latitude=country_center_y,
                                       longitude=country_center_x)
        if country_center_y < 0:  # if point lies south of equator
            self.crs_dict = {'proj': 'utm', 'zone': specific_utm[2], 'south': True}
        else:  # if point lies north of equator
            self.crs_dict = {'proj': 'utm', 'zone': specific_utm[2]}

To Consider

I do not have a feeling for how much off the calculation of the centroids are close to the poles. Mabye it is negligible if the centroid is calculated on a country or 10×10 degree tile base.

Best Answer

With recent (>= 0.9) versions of geopandas, GeoDataFrame has an estimate_utm method that should do just that, automatically.

Related Solutions

[GIS] How to Display the World Map in Projected Coordinate Systems

Edit:
The problem is that ALL projected systems are made to...

Give accurate measurements.
or
Give a pretty picture.

An accurate measurement projection is for a small"ish" area. Not the entire world.
All this being said. Arcgis will allow (other software do this too) multiple projections into one or different projections.
So you just decide what projection "looks" the best for your map, Set the document projection (see previous answer from below) to that. and then load all your data.
The only reqirement is that each data source would need to have a projection defined. (shape files use a *.prj file, tif images use a *.tfw file).

These helper files tell the software what the coordinates relate to. On-the Fly projection does the rest of the work.

I suspect that because of the question you have one or more datsesets that are not overlaying correctly.
If this is true you should look at the prj file to determine if: first there is one; and second if it is defined as the correct projection.
You can search gis.stackexchange for other questions related. LIKE: "projections"
Or just look at the related to the right of this question.

If you look at the TOC (table of contents). You should have a data frame named layers (that is just a default name).
Right click on that and select properties.
Go to the coordinate systems tab and select the UTM zone you want.
CRS
Unless you are actually asking about how to show the utm grid (which is in the grid tab [same dialog])
this should do what you want.

[GIS] Clipping Rasters from Mosaic Dataset in Original Coordinate System

Set the Output Coordinate System environment setting when clipping the raster using Extract by Mask with an input shape file. If you are implementing this in python you should be able to grab the coordinate system from the shape file using the describe method.

Projecting/converting WGS84 World Mercator to WGS84 UTM zone 14N does not require a parameter transformation as they are based on the same geographic coordinate system, if they were dissimilar then this would not work.

Here is a script that will clip a raster mosaic to a shapefile and match the output coordinate system to the shapefile:

import arcpy, sys

InRasterMosaic = sys.argv[1]
InShapeFile = sys.argv[2]
OutRaster = sys.argv[3]

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddMessage("Unable to get spatial analyst extension")
    sys.exit(0)


desc = arcpy.Describe(InShapeFile)
SR = desc.spatialReference
arcpy.env.outputCoordinateSystem = SR

ExRas = arcpy.sa.ExtractByMask(InRasterMosaic,InShapeFile)
ExRas.save(OutRaster)
arcpy.CheckInExtension("Spatial")

I have tested this and it works! Have a read of Adding a Script Tool, follow those instructions and you can put this into a model.

If the input is in WGS (DD) as stated:

import arcpy, sys, os

InRasterMosaic = sys.argv[1]
InShapeFile = sys.argv[2]
OutRaster = sys.argv[3]

if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage("Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddMessage("Unable to get spatial analyst extension")
    sys.exit(0)

# from the geographic extent, get the centre
# of the extent, only the X is required
desc = arcpy.Describe(InShapeFile)
Ext = desc.Extent
CentX = ( Ext.XMin + Ext.XMax ) / 2

# I'm not familiar with North American projections
# so you'll have to put the numbers in for each zone
if CentX < Zone14max:
    ProjSR = arcpy.SpatialReference(32614)
elif CentX < Zone15max:
    ProjSR = arcpy.SpatialReference(32615)
# continue for the required spatial references
else:
    arcpy.AddMessage("Unable to establish coordinate system")
    arcpy.CheckInExtension("Spatial")
    sys.exit(-1)

arcpy.env.outputCoordinateSystem = ProjSR

ExRas = arcpy.sa.ExtractByMask(InRasterMosaic,InShapeFile)
ExRas.save(OutRaster)

arcpy.CheckInExtension("Spatial")

This will need you to make some edits, but will basically do what you're after. Using the centroid of the extent work out which zone it's likely to be in and then use that for the output coordinate system. The numbers are the SRID of the output coordinate system like this one http://spatialreference.org/ref/epsg/32614/

Best Answer

Related Solutions

[GIS] How to Display the World Map in Projected Coordinate Systems

[GIS] Clipping Rasters from Mosaic Dataset in Original Coordinate System

Related Question