Clipping in Antarctica #278

carocamargo · 2021-03-23T18:47:01Z

carocamargo
Mar 23, 2021

Hi,
I have a shapefile of drainage basins from Antarctica, and a global ocean-land mask.
I am trying to select the regions in my mask that correspond to the polygons in the shapefile, but without success so far.
Any suggestions?
Ideally, I would like to set different numbers to the different geometries of the shapefile.

This is what I got so far:

import geopandas as gpd
import numpy as np
import rioxarray
import xarray as xr
#%%
path='/usr/data/'
ds=xr.open_dataset(path+'LAND_MASK_CRI-JPL_180x360_conservative.nc')
ifile='ANT_Basins_IMBIE2_v1.6'
shp=gpd.read_file(path+'ANT_basins_shp/'+ifile+'.shp')
print(shp.type)
print(shp.crs)
print(shp.bounds)
#%%
bounds=np.array(shp.bounds)
xmin=bounds[1:19,0]
ymin=bounds[1:19,1]
xmax=bounds[1:19,2]
ymax=bounds[1:19,3]
#%%
shp.plot()
#%%
da = ds.mask
# Convert lon to -180 - 180
da = da.assign_coords(lon=(((da.lon + 180) % 360) - 180))
rda = ds.rio.set_spatial_dims('lon', 'lat')
rda = rda.rio.write_crs("epsg:4326")
clipped = rda.rio.clip(shp.geometry, "epsg:4326")
clipped.mask.plot()

Thank you

ANT_basins_shp.zip
LAND_MASK_CRI-JPL_180x360_conservative.nc.zip

Answered by carocamargo

Mar 26, 2021

Finally got it working.. I think the issue was with the shapefile data itself. The coordinates of the shapefile were in '3031' projection, and just switching the proj=3031 in clipping wasn't working, so I had to reproject the shapefiles before hand.
not sure if it's the best solution, but it's working.

import pandas as pd
# import geopandas as gpd
import numpy as np
import rioxarray
import xarray as xr
from cartopy import crs as ccrs#, feature as cfeature
from matplotlib import pyplot as plt
import matplotlib.path as mpath

from pyproj import Transformer
import shapefile 
#%% Open 1 degree mask
path='/Users/ANT_basins_shp/'
ifile='LAND_MASK_CRI-JPL_180x360_conservative.nc'
ds=xr.open_data…

View full answer

snowman2 · 2021-03-23T20:20:27Z

snowman2
Mar 23, 2021
Maintainer

#230 may be helpful

1 reply

carocamargo Mar 24, 2021
Author

Thanks, if shone some light... but still not working wit the shapefiles itself...
If I use box from shapely, it works well:

xda = xr.DataArray(data=ds.mask, 
                       dims=["y", "x"], 
                       coords=dict(x=(["x"], ds.lon),
                                   y=(["y"], ds.lat)))
# transform long form 0->360 to -180->180
xda = xda.assign_coords(x=(((xda.x + 180) % 360) - 180))

proj="EPSG:3031" #EPSG:3031: WGS 84 / Antarctic Polar Stereographic
# proj="EPSG:4326"

xda=xda.rio.write_crs(proj)

geometries = [box(-180, -90, 180, -70)]
clipped = xda.rio.clip(geometries, crs=proj)

But then, when i try the geometry form the shapefile, it says there is no data...

xda = xr.DataArray(data=ds.mask, 
                       dims=["y", "x"], 
                       coords=dict(x=(["x"], ds.lon),
                                   y=(["y"], ds.lat)))
xda = xda.assign_coords(x=(((xda.x + 180) % 360) - 180))

proj="EPSG:3031" #EPSG:3031: WGS 84 / Antarctic Polar Stereographic
# proj="EPSG:4326"

xda=xda.rio.write_crs(proj)

geometries = [shp.geometry[13]]
clipped = xda.rio.clip(geometries, crs=proj)

This is the error I get:

And I tried both with the normal projection (4326) and with the Polar one (3031).

My shapefile is in the polar projection. If I do:

print(shp.crs)

{'init': 'epsg:3031'}

and:

shp.plot()

I tried converting the latitude and longitude of the mask to polar projection, but it still didn't work as well...
Any ideas/suggestions?

carocamargo · 2021-03-26T15:09:21Z

carocamargo
Mar 26, 2021
Author

Finally got it working.. I think the issue was with the shapefile data itself. The coordinates of the shapefile were in '3031' projection, and just switching the proj=3031 in clipping wasn't working, so I had to reproject the shapefiles before hand.
not sure if it's the best solution, but it's working.

import pandas as pd
# import geopandas as gpd
import numpy as np
import rioxarray
import xarray as xr
from cartopy import crs as ccrs#, feature as cfeature
from matplotlib import pyplot as plt
import matplotlib.path as mpath

from pyproj import Transformer
import shapefile 
#%% Open 1 degree mask
path='/Users/ANT_basins_shp/'
ifile='LAND_MASK_CRI-JPL_180x360_conservative.nc'
ds=xr.open_dataset(path+ifile)
ds

xda = xr.DataArray(data=ds.mask, 
                       dims=["y", "x"], 
                       coords=dict(x=(["x"], ds.lon),
                                   y=(["y"], ds.lat)))
# Convert longitudes to -180 to 180
xda = xda.assign_coords(x=(((xda.x + 180) % 360) - 180))
# # Order longitudes 
xda = xda.sortby('x')
xda.plot()
#% %
# Set the projection
proj="EPSG:3031" #EPSG:3031: WGS 84 / Antarctic Polar Stereographic
# proj="EPSG:4326"
xda=xda.rio.write_crs(proj)
# xda.rio.set_crs("EPSG:3031")
#%% open shapefile
ifile='ANT_Basins_IMBIE2_v1.6'
sf = shapefile.Reader(path+ifile+'.shp')
shapes = sf.shapes()
#%% Transformation function
inProj=3031#EPSG:3031: WGS 84 / Antarctic Polar Stereographic
outProj =4326
transformer = Transformer.from_crs(inProj, outProj)
for i in range(len(shapes)):
    #shape[0] represents the 1st polygon of the shapefile
    # for point in shapes[0].points:
    #     print (point)
    # Get coordinates of each polygon:
    coords=[point for point in shapes[i].points]
    #Transform the coordinates from 3031 to 4326
    points=[pt for pt in transformer.itransform(coords)]
    # lat,lon = transformer.transform((lat,lon))
    x=np.array([p[1] for p in points])
    y=np.array([p[0] for p in points])
    df2=pd.DataFrame(x,columns=['lon'])
    df2['lat']=y
    coords = df2.values.tolist()
    # Build the geometries object
    geometries = [{'type': 'Polygon','coordinates': [coords]}]
     clipped = xda.rio.clip(geometries)

The final result, after correcting a bit to include the center of the pole:

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Clipping in Antarctica #278

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 2 comments 1 reply

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

Clipping in Antarctica #278

Uh oh!

Uh oh!

carocamargo Mar 23, 2021

Replies: 2 comments · 1 reply

Uh oh!

snowman2 Mar 23, 2021 Maintainer

Uh oh!

carocamargo Mar 24, 2021 Author

Uh oh!

carocamargo Mar 26, 2021 Author

carocamargo
Mar 23, 2021

Replies: 2 comments 1 reply

snowman2
Mar 23, 2021
Maintainer

carocamargo Mar 24, 2021
Author

carocamargo
Mar 26, 2021
Author