texts: How to extract CRS ellipsiod info in correct format for pyproj inverse

texts:
I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe (gdf), but my datasets are in different locations around the world. I need to use a local UTM EPSG to get accurate azimuths and distances (e.g., discussion here; this is common knowledge).
For a given EPSG, how can I systematically extract the g = pyproj.Geod(ellps='X') info from the CRS in the right format for X?
Below is my best attempt using myellipsoid = CRS.from_user_input(myepsg).ellipsoid, but it's in the wrong format. In this example, "GRS 1980" needs to be "GRS80"
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS

myid = [1, 1, 1]
myorder = [1, 2, 3]
lat = [5174925.07851924, 5174890.26832387, 5174855.45812849]
long = [1521631.6994673, 1521667.11033893, 1521702.52121056]
myepsg = 2193

df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=myepsg)
myellipsoid = CRS.from_user_input(myepsg).ellipsoid

print(myellipsoid)
print(gdf_pt.crs)
display(gdf_pt)

ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);

g = pyproj.Geod(ellps=myellipsoid)
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = g.inv(gdf_pt.long[i], gdf_pt.lat[i], gdf_pt.long[i-1], gdf_pt.lat[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist'] = myinfo[2]
    gdf_pt.loc[i, 'bearing'] = max(myinfo[1], myinfo[0])

display(gdf_pt)

Using: Windows 10; conda 4.8.2; Python 3.8.3; shapely 1.7.0 py38hbf43935_3 conda-forge; pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge


texts: How to change or assign ArcMap layer symbology type from python?

texts:
How do I SET an ArcMap layer's symbology type using python?

Layer()SymbologyType property is read only and UpdateLayer() requires a pre-existing layer source to copy symbology from. Copying from another source doesn't work well because the properties to set are dynamic. For example derive break values for classified raster from the input source and an algorithm.

All I want to do is assign the type, e.g.

lyr.symbologyType = "RASTER_CLASSIFIED"
#or
lyr.update(symbologyType="UNIQUE_VALUES")


texts: Making GeoServer run as Windows Service after being installed?

texts:
While installing GeoServer, on my windows machine, I selected the 'Run manually' as the option.

Now I wish to keep it running all the time, and for this, I would have to run it as a Windows service.

I know that I can manually make a new service and run it from there, but is there a better way to make it run as a Service, without reinstalling it?


texts: How to build a relational database system with 4 tables in arcgis desktop?

texts:
I have been given data from an working online application which consists of a relational model. My aim is to recreate the model using Arcgis Desktop.

I have the table "elements" containing Point geometry with the attribute element_id as the primary key. In the table Aspects, Types and media. The tables "Aspects" and "Media" have a 1:1 relationship with "elements" and media has a 1:m relationship.

It is possible to build these multiple relationships in ArcGIS Desktop?


texts: Shaded relief from DEM in GRASS and QGIS

texts:
I have been trying to get a good looking shaded relief model out of GRASS and Qgis for a day or so now and just can't seem to get it.  The biggest problem I'm having is the valleys look like mountains, but I just can't seem to get the right combination in order to show mountains and valleys look right.

Bellow is a some pictures of the problems:

DEM (Brown is high and Red is low, red is only used for a visual check)

Shaded Model 1 with DEM Sun position 180 degress from North, sun height 80

Problem: Hills look good but the river (in red) still looks like a mountain even though it's the lowest point on the map.

Shaded Model 2 with DEM Sun position north west (315 degrees), sun height 60

Problem: Again the hills look good but the river (in red) still looks like a mountain even though it's the lowest point on the map.

I have tried all different combinations of angles and sun heights and I just can't seem to get the river to look like a river, it always looks like a hill.

Can anyone point me in the right direction on what I might be doing wrong?

EDIT: So I took the contour data home and ran the process from the start and used the defaults for r.shaded.relief and got this:

This looks much better to me, not the greatest colours but we can change that.  I think my work install of grass must be busted because even the defaults looked bad.


texts:
There are some similar question here, but nothing worked for me. So I have a database that has the ending .mdb which to my knowledge is some kind of Microsoft access database. But I also tried it with a test-database that has the standard .accdb-ending and did not manage to make it work.
I did the following: I clicked on Layer --> addd Layer --> add Vector Layer. Then I select database and type=ODBC and opening the connection dialog this appears:

And I'm not really sure what to enter here. I tried several thigs and none of them did work for me. So basically I guessed that I can name it as I want to name it. And the most important part is the Database-field. However, when I put the path to the database in there I always get the error:
Connection failed - Check settings and try again.

Extended error information:
Unable to initialize ODBC connection to DSN for C


texts: I am a B.Tech graduate(2013) in Geo Informatics and Geo Technology.Presently, i am working i an company as Junoir Executive in Gis . I also know some  programming languages like "V.B.net" and C#.My doubt is how can i improve my Career strength in this field so that i can earn handful salary ?


texts: Environmental Analytical Lab Data to Geodatabase

texts:
What is the best way to get groundwater analytical data into a GIS geodatabase?

The data is currently in Excel format, but there are options from the lab for other file extensions dbf.

What is the best way to format or import analytical data from a groundwater monitoring well event?


texts: How to zoom out in Openlayer in the particular extent?

texts:
I am creating the interaction for zoom in and zoom out using OpenLayers. I am able to achieve zoom in but I am not able to figure out how can I create zoom out to some extent. Here is my code for zoom in
$('#zoomin').click(function () {
    removeInteraction();
    $('#map').css("cursor", "zoom-in");
    map.addInteraction(dragBoxInteraction);
    dragBoxInteraction.on('boxend', function () {
        var extent = dragBoxInteraction.getGeometry().getExtent();
        map.getView().fit(extent, map.getSize());
    })
})

For zoom out I have written something like this
$('#zoomout').click(function () {
    removeInteraction();
    map.getView().animate({
        zoom: map.getView().getZoom() - 1,
        duration: 250
    })
})

Using this way I am able to draw a rectangle box and I am getting the extent and I am able to zoom in that extent.
But Using the same feature how will I zoom out?


texts:
I have a water feature layer which includes rivers, lakes and streams. I want to create, essentially, a buffer of these features equal to that of a 50cm increase in elevation from its current boundary. I have been toiling with this one for a few hours now, but can't theoretically get through the problem.

The result would be a buffer that is "aware" of the elevation in its determination of the buffer edge.

Any way I think about this, I keep coming back to a looping function addressing each point along the boundary of each water feature. This would be computationally intensive and I don't have all the theory down for that method either.

What do you think?

ArcGIS 10.1, all license levels with LiDAR DEM.

The water is, itself, a constant elevation, but as the land rises, so too does the channel. If this isn't clear, think of a mountain stream flowing to the valley.

Have successfully made a layer from suggestions below by Michael.

The result is not ideal (used a 1m elevational increase), thus looking to see if anyone can improve upon this method.


texts: Generating new data in curve (polyline) by interpolation of existing data in this curve

texts:
I am trying to add altitude data on the vertex of a polyline in QGIS.
I have extracted points from the vector layer of the polyline to update the Z-value.
I have altitude for some points, new points have been added by a refinement of the initial curve.
New points have a NULL value for altitude ("elev_ft").

How can I add altitude for these new points by interpolation on the altitudes already known?
Interpolation should be based on the distance of the path between the points.


texts: Plotting Shapely Multipolygon using Matplotlib

texts:
I have few polygons or boxes. Ones intersect each other but some are isolated. For example I have free figures:
import shapely.geometry as sg
import shapely.ops as so
import matplotlib.pyplot as plt

r1 = sg.Polygon([(0,0),(0,1),(1,1),(1,0),(0,0)])
r2 = sg.box(0.5,0.5,1.5,1.5)
r3 = sg.box(4,4,5,5)

First two are intersected and r3 is on some distance from them. I unite them via cascaded_union:
new_shape = so.cascaded_union([r1, r2, r3])

Then I try to plot it (one united figure of r1 and r2 and  one distanced box r3)
xs, ys = new_shape.exterior.xy
fig, axs = plt.subplots()
axs.fill(xs, ys, alpha=0.5, fc='r', ec='none')
plt.show()

Instead of plot I receive an AttributeError: 'MultiPolygon' object has no attribute 'exterior'.
Is there a pythonic way to display a multipolygon  or to iterate through it and draw its parts?


texts: Is postgres 9.1.9 supported by ArcGis 10.1 Service Pack 1?

texts:
I want to install Postgres and Postgis.

Now we use ArcGis 10.1 So I install Postgres 9.0.5 But there is not exist Postgis for this version. I try install postgis 2.0.3  But I fail(it is only 32 bit application ) http://download.osgeo.org/postgis/windows/pg90/

So I'd like to install postgres 9.1.3 But I find only 9.1.9 version :(

I'd like to install Postgis because we will use geoserver.


texts: Referencing Bookmarks in Alphabetical Order in Python script tool?

texts:
I have a custom script tool that is supposed to set up a map layout, by adjusting the scale bar, north arrow, title, and legend. It is then supposed to reference a bookmark and export the layout as a PDF.

The problem is that when I use this tool in ModelBuilder it references the wrong bookmarks

I'm using this script tool in ModelBuilder for ArcGIS Pro. This model iterates through Urgent Care centers and its corresponding data one at a time (in alphabetical order by center name). I have bookmarks for every Urgent Care center, so my thought is that, if I can find a way for the model to reference the bookmarks in alphabetical order, then they will correspond with the correct Urgent Care center.

Could this be achieved via a table or list Index, or some kind of loop?
Or is there a better way to go about this whole process?

Here is the code from the script tool I am using.

import arcpy

# input layer

lyr = arcpy.GetParameterAsText(0)

# input name of layout

p = arcpy.mp.ArcGISProject("CURRENT")
lyt = p.listLayouts("Layout_King")[0]

# Reposition the scale bar

scaleBar = lyt.listElements("MAPSURROUND_ELEMENT", "Scale Bar")[0]
mf = scaleBar.mapFrame
scaleBar.elementPositionX = mf.elementPositionX + 0.0
scaleBar.elementPositionY = mf.elementPositionY - 0.5

# Reposition the north arrow

northArrow = lyt.listElements("MAPSURROUND_ELEMENT", "North Arrow")[0]
mf = northArrow.mapFrame
northArrow.elementPositionX = mf.elementPositionX + 8.8
northArrow.elementPositionY = mf.elementPositionY + 0.7

# Align the title with the center of the map frame

title = lyt.listElements("TEXT_ELEMENT","Name of Map Text")[0]
mf = lyt.listElements('MAPFRAME_ELEMENT',"Map Frame")[0]
title.elementPositionX = mf.elementPositionX + (mf.elementWidth / 3.7)
title.elementPositionY = mf.elementPositionY + (mf.elementHeight / 0.98)

# Reposition the Legend and fix legend title

legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
legend.title = "Legend"
legend.elementPositionX = mf.elementPositionX + 7.7
legend.elementPositionY = mf.elementPositionY + 7.15

# setting layout to bookmark

aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout

lyt = aprx.listLayouts("Layout_King")[0]
mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark

bkmks = mf.map.listBookmarks("*")
for bkmk in bkmks:
  mf.zoomToBookmark(bkmk)
lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")

I tried the bkmks.sort() earlier, but it gives me this error:

Traceback (most recent call last):
  File "<string>", line 10, in <module>
TypeError: unorderable types: Bookmark() < Bookmark()

I also tried to print the bookmark list using:

bkmks = mf.map.listBookmarks("*") bkmks.sort() print (bkmks)

But it gives me this error:

 File "<string>", line 1
    bkmks = mf.map.listBookmarks("*") bkmks.sort() print (bkmks)
                                          ^
SyntaxError: invalid syntax


texts:
Is there a source (free or paid) for polygons of border coordinates for spanish postcodes?
Google show this borders in google maps, but doesn't provide them.


texts: Web-based user editable routing application

texts:
I am totally new to web-based GIS (and open source GIS) and am now working on a web-based travel planner using Postgis and Pgrouting with multicriteria costs. I want the end users (the traveler) to be able to decide the cost for different road links.

So I need a web-map service which has vector road links that will be possible to click and edit for logged in users. The edits will then be saved to the postgis database, and the average cost be updated.

Any recommendations for solutions to this? Is Geoserver a good way to go? Mapfish? Is this a tricky thing or easily done (I am not a programmer or system developer, but willing to learn, though)?

I´d really appreciate any tips and help with this


texts: Different number of grid cells between years when using a grid layer to look at area of land cover change over time

texts:
My goal is to calculate the area of mangrove cover (vector layer) per grid cell (vector layer) over time (yearly), to determine land cover change.

I have tried two ways so far.
Initially I wanted to Intersect the grid and polygon layers then dissolve by ID and then calculate the area. This has two problems 1) it does not include any grid cells with 0 mangrove polygon intercept, and 2) there are duplicated grid ID values where two polygons intercept the same grid.

I also tried Union. This includes the grid cells with 0 mangrove polygon, but still has the issue with duplication.

Any ideas on how I can avoid this duplication so that I have the same number of grid cells per year so that I can compare between years?


texts: Why does the flow direction tool in ArcGIS desktop give wrong output (unique values from 1-255)?

texts: I ran the flow direction in ArcGIS, and received unique values ranging from 1-255. I know that these values are wrong, and should be along the lines of 1, 2, 4, 8, 16, etc multiples. The DEM is projected. What could cause this?


texts: PyQGIS: How can I load a specific layer from a file containing multiple

texts:
I have a vector file that contains multiple layers, if I open it in QGIS I get the "Select vector layers to add dialog" and I can choose one and it loads fine.

However in a .py script if I load it via

layer = QgsVectorLayer(filename, "mylayer", "ogr")

I just get one of the layers (layer 0)

How do I specify which layer I want to load?
Or if it IS actually loading all of them, how do I change to the one I want?

My full script

import os
import sys
import qgis
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import ogr

qgs = QgsApplication([], False)
QgsApplication.setPrefixPath( r"C:\Program Files\QGIS 2.18\apps\qgis-ltr", True )
QgsApplication.initQgis()

filename = "C:/datafiles/sourcefile"
if os.path.exists(filename):
    print('found it')

layer = QgsVectorLayer(filename, "mylayer", "ogr")
if layer.isValid():
    print("Layer was loaded successfully!")
else:
    print("Layer failed to load!")
    sys.exit()

print (layer.type())
print (layer.name())
print (layer.id())

print('finished')
qgs.exitQgis()


texts: ImportError: No module named qgis.core on Windows 10

texts:
I am trying to execute the following code on my Windows 10 with QGIS 2.18.20.

from qgis.core import *
from qgis.utils import *
import os, sys

QgsApplication.setPrefixPath("C:/Program Files/QGIS 2.18/apps/qgis-ltr", True)
app = QApplication([], True)
QgsApplication.initQgis()

sys.path.append('C:/Program Files/QGIS 2.18/apps/qgis-ltr')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

path_to_layers = "path"

gebaude = path_to_layers + "layer1.shp"
duek = path_to_layers + "layer2.shp"
ergebnisse = path_to_layers + "result.shp"

general.runalg("qgis:joinattributesbylocation", gebaude, duek, ['intersects','overlaps','within'], 0, 0, 'sum', 1, ergebnisse)

print("Finished with no issues")

However, I am encountering the error

Traceback (most recent call last):
  File "path_to_the_python_code\Gebaeude.py", line 1, in <module>
    from qgis.core import *
ImportError: No module named qgis.core

Have already seen similar threads on this topic but can not really understand how and where I should adjust the full path to the Python executable under the QGIS installation.

set PYTHONPATH="C:/Program Files/QGIS 2.18/apps/Python27"
set PATH="C:/Program Files/QGIS 2.18/bin"

I tried the above set-ups in OSGeo4W Shell through the bat file but no success yet. What I am doing wrong?

References:

Fixing “ImportError: No module named qgis.core”?
Problem with import qgis.core when writing a stand-alone PyQGIS script


texts: How to save data from QGIS 1.8 to PostGIS without shapefile intermediate?

texts:
Is it possible to save data directly into an existing PostGIS database from QGIS without having to first save it as a shapefile? I know about the SPIT plugin but that requires a shapefile. I want to know if is there a way to skip the save-to-shapefile-use-SPIT/shapefile-loader-to-import-to-Postgis process after an analysis in QGIS.

Edit:
Example: Load road layer from PostGIS to QGIS. Calculate shortest path using road graph plugin and export as temporary (temp) layer. Right click on temp layer, select save as... the save as options include the general vector type files (shapefile,GeoJSON) and at the bottom of the list is an option to save as a spatiallite database. This, unfortunately, only creates a new spatiallite database with only the one shortest path road layer (no postgis option). What I would like to know is, is it possible to save a temp file directly into an existing PostGIS database?


texts: Shapefile saving to unknown format instead of geojson

texts:
I have a shapefile and would like to export a single record (polygon), which I've selected, to GeoJSON.

I have never had a problem with this before.  I've checked the CRS - it is defined as ESPG: 4326.  I choose Save As... > Tick "only selected features" but then I am given an unknown table instead of a GeoJSON file.

Why might this be happening?

Here is a link to the shapefile in question: link


texts: Excluding field on table from being copied by ArcPy?

texts:
I have a query "layer":

>>> r = arcpy.MakeQueryLayer_management(r'C:\path\to\connection\file\connection.sde', 'TEST', 'SELECT T.ID AS OBJECTID, T.* FROM MY_TABLE T WHERE T.SOME_COLUMN = 1', oid_fields=['OBJECTID'])
>>> r
<Result 'TEST'>
>>> r[0], type(r[0])
(u'TEST', <type 'unicode'>)
>>> table = arcpy.mapping.TableView(r[0])
>>> len(arcpy.ListFields(table))
22

Since the database table does not have a spatial column, MakeQueryLayer_management returns a string that can then be turned into TableView.

I want to copy the features into a file geodatabase, but I want to exclude SOME_COLUMN from the output. Since the documentation for MakeTableView says

      A subset of fields can be made unavailable in the new layer by using the Field Info control's visible property. ... You cannot use the hidden fields in a workflow if the newly created layer is input to a subsequent process or tool. If the output is saved to disk, only the fields listed as visible will appear in the new data,

it seemed like the perfect fit. But I'm getting strange results, and the field is showing up in my final output anyway:

>>> field_changes = arcpy.FieldInfo()
>>> field_changes.addField('SOME_COLUMN', 'SOME_COLUMN', 'HIDDEN', 'NONE')
>>> new_table = arcpy.MakeTableView_management(table, 'NEWTEST', field_info=field_changes)
>>> new_table
<Result 'NEWTEST'>
>>> len(arcpy.ListFields(new_table))
22

So this is where it starts getting weird. It looks like all the columns are still on the new table. SOME_COLUMN is the last column on my table, so let's take a look:

>>> some_column = arcpy.ListFields(new_table)[-1]
>>> some_column
<Field object at 0x12a7c1f0[0x29e02a8]>
>>> some_column.name
u''
>>> some_column.type
u'Integer'
>>> some_column.length
4

...?  (I have confirmed that all the other fields are account for, so this Field object does correspond to the field I tried to hide.) Apparently, the field is still there, just with an empty string name. When I copy it, the field gets its name back and shows up in the final output:

>>> arcpy.CopyRows_management(new_table, 'C:/path/to/mygdb.gdb/TESTOUT')
<Result 'C:\\path\\to\\mygdb.gdb\\TESTOUT'>

The funny thing is that I can make it work with actual Layers that have a spatial column. It still does the weirdness with the field having an empty name, but the field doesn't show up in the actual output:

>>> layer = arcpy.MakeQueryLayer_management(r'C:\path\to\connection\file\connection.sde', 'TESTLAYER1', 'SELECT T.ID AS OBJECTID, T.* FROM MY_SPATIAL_TABLE T WHERE T.SOME_COLUMN = 1', ['OBJECTID'])[0]
>>> layer
<map layer u'TESTLAYER1'>
>>> len(arcpy.ListFields(layer))
17
>>> new_layer = arcpy.MakeFeatureLayer_management(layer, 'TESTLAYER2', field_info=field_changes)[0]
>>> new_layer
<map layer u'TESTLAYER2'>
>>> len(arcpy.ListFields(new_layer))
17
>>> x = arcpy.ListFields(new_layer)[-1]
>>> x
<Field object at 0x12e7ebf0[0x29e07b8]>
>>> x.name
u''
>>> arcpy.CopyFeatures_management(new_layer, r'C:\path\to\mygdb.gdb\NEWLAYER')
<Result 'C:\\path\\to\\mygdb.gdb\\NEWLAYER'>

What am I doing wrong?

Extra Info:

I tried MakeQueryTable since it's also supposed to be able to give me a subset of fields, but that gave me the error ERROR 000055: Cannot create a Query Table for this workspace.

I could just create the table and then delete the field after creating it, but that would force more of my code to be aware of the filtering column (it's name, existence, how it's handled). I'd prefer to centralize all that handling into the code that actually creates the query layer.


texts: How to convert gml:pos to latitude and longitude?

texts:
I'm trying to figure out the latitude and longitude of rail stations that are contained in this WMS feature list. I'm completely new to the GIS field, so please forgive me if I use the wrong terms below. I've done quite a bit of reverse engineering and fiddling to get this far.

The first station has the following location:

<gml:Point srsName="urn:ogc:def:crs:EPSG::27700" srsDimension="2">
  <gml:pos>450077 356089</gml:pos>
</gml:Point>

After some research, I've concluded this position is relative to offsets and I've used the EPSG site to look up code 27700. This shows me:

  Latitude of natural origin: 49° N
  Longitude of natural origin: 2° W
  Scale factor at natural origin: 0.9996012717
  False easting: 400,000 (metre)
  False northing: 100,000 (metre)

I'm guessing the values in my station location have had 400,000 and 100,000 respectively added to them. So if I subtract those values I would get:

50077, 256089

But I'm not sure what these values represent. They are minuscule, so they can't be metres from the latitude/logitude point (0,0).

Can anyone point out where I've gone wrong? My goal is to eventually convert this gml:pos value into latitude and longitude.


texts: How to emulate OpenLayers button click in C#

texts:
I'm trying to call a click of a user using C # with the following code ( " C # :

  GeckoWebBrowser.Navigate ( javascript : void (
  Document.getElementById ( ' idOfElement ' ) .click ( ) ) " ) ;

But this element is not triggered in Open API Layer .

How can I trigger an action using this code snippet ?

What is the function calling this action ?


texts: Editing multiple feature classes contained in feature dataset using ArcMap?

texts:
I have two feature classes contained inside a file geodatabase feature dataset.
One is point and one is polygon.
When I turn on editing mode for one of them, the other is editable as well.

Assuming that is the correct behavior by ArcMap, is that one of the advantages of using a feature dataset?


texts: How to create themed layer based on a list of geometries and specific risk values in Sharpmap

texts:
I am now using SharpMap for creating a network-based tool. By using SharpMap, for constructing network, the current functions in SharpMap enables the tool to capture coordinate on background layer and then creates geometrical points and polylines. Following this step I could create  either point-based, or polyline-based vector layers.

Now I need to display each geometry into different colours on same vector layer based on their "risk index": a value to assess the gradient colours needed to use. From SharpMap discussion I understand we could not add data attribute to the layers in the 1.0 version, and only Data Column in vector layer can be used to render the colours of elements.

My difficulty is: since I created the vector layer by myself and not imported from external data source, and therefore there is no attribute table attached to each layer, it means, no attribute table in the layer. Therefore, I am wondering are there any methods available that I can render my layers not based on attribute table?


texts: Date of reference in an feature collection

texts:
I am pretty new on GEE and idk how easy is this task. I want to find daily precipitation for all Brazilian municipalities in 2006. I applied a map function with reduce regions and end up with a feature collection of municipalities. But I really need the date of reference for each element (day and month). I tried to use set function but idk how to guarantee each information would have the correct date from the image that was calculated. My ending goal is a table to export to google drive.

Here's my code:

var mun = ee.FeatureCollection("users/rodrigoccuri/BRMUE250GC_SIR");

// Importing Chirps dataset:
var chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY');

// Filtering year 2006:
var chirps06 = chirps.filterDate('2006-01-01', '2006-12-31');

// Creating a function to reduce regions for each image by pixel mean for mun:
var reduce = function(image) {
  return image.reduceRegions(mun, ee.Reducer.mean(), 1000)
}

// Appling reducer for each image in the collection:
var precmun06 = chirps06.map(reduce).flatten();

print(precmun06.first());


texts: Free GIS workshops, tutorials, and applied learning material

texts:
This Q&A lists free as in $0 workshops in GIS and related fields that have their material available to view or download online. The list is not limited to FOSS GIS, as GISers usually use a combination of open and closed source tools to accomplish their tasks. Some of these tools are easier than others and some are better documented. This list will be a great benefit to the community especially for new comers to the field.
General GIS

GIS Project Video Tutorial on Acquiring, Analyzing, and Mapping US Census Data in QGIS from @A.S
GIS Tutor: Beginner / Intermediate Level GIS from @radek
BostonGIS from @simo
Lectures on GIS for the Social Sciences from @ubernatural

Open-source GIS

OpenGEO Education Center from @radek
QGIS Videos (faunitalia) from @simo
GIS SE Question: Geoserver Tutorials from @com
GIS SE Question: Open Source Training Materials from @MarkIreland
QGIS for Newbies from @IanAllan

ESRI Products

GIS SE Question: Best place for (structured) ArcGIS tutorials from @robintw
ArcGIS Automation and Programming from @Bethany

Transportation
These links include either theoretical or applied transportation knowledge in transportation planning

GIS primer for Transportation
Transit Capacity and Quality of Service Manual
UrbanSim: FOSS Urban development, socio economic, and land use planning package
MOVES(Motor Vehicle Emission Simulator) Workshop and technical background
Transit GIS tutorials
Remote sensing in transportation workshop
Florida's CBT planning model explained

Disaster Management

Geo-information Technology for Crisis Management

Databases

Spatial Database Course Material from @radek

Geo Statistics

Learn R Lectures and Classes
Geodatabase Mining Course from @radek
A Practical Guide to Geostatistical Mapping


texts: I have a shapefile where in I want to select features based on their feature ids

texts:
I have a list of feature ids and I am using layer.SetAttributeFilter for that but the code neither selects those features nor does it raise an error. I took help from these questions:

OGR Filter Vector Layer by a List of Feature ID's

OGR: How to save layer from AttributeFilter to a shapefile?

Here is my code:

inlyr=inds.GetLayer()
m =[]
for a in range (0, len(new_coords1)):
    m.append(new_coords1[a])
inlyr.SetAttributeFilter("POINTID IN {}".format(tuple(m)))
drv = ogr.GetDriverByName( 'ESRI Shapefile' )
outds = drv.CreateDataSource( "test2.shp" )
outlyr = outds.CopyLayer(inlyr,'test2')
del inlyr,inds,outlyr,outds

where POINTID is the name of the id field. The output creates a new shapefile with all the features and not the selected features. I also tried executing the SetAttributeFilter separately which gives an output of 5. I don't know what that output means.


texts: Moving "last name" from last to first in Attribute Table of ArcGIS Desktop?

texts:
How do you get e.g. "John D. Doe" to "Doe, John D." in the AT using VB Scripting without going individually going through each line?  I am using ArcGIS 9.3.

I should note that not all the names are the same; they are different.  Is it possible to select those in the column with different names and just reverse the "first, last name" to "last name, first"?


texts: Show labels for point features outside map canvas

texts: I'm using a pyqgis script to save multiple images of a map and then stitch the images together to create a larger image. One of the layers is a point layer containing all the labels. I noticed that if a point is outside the current map canvas extent, the label for that point is not shown even if the actual label text falls within the extent of the map canvas. This results in labels being cut off at the edges of the saved images. Is there a way to force labels to be shown even if the feature falls outside the map canvas extent?


texts: What are the rules for releasing ArcObjects from memory in .NET?

texts: Do I need to release each object?  Is some memory management handled for me?


texts: Find polygon boundary dynamically based on a point (lat/long) having type farmland using OSM

texts:
There are many ways to retrieve the polygon information when you know that the information is already present somewhere (data locally or in database containing extracted info about polygon). But this question is to figure out, how do we retrieve a polygon info when you don't have any extracted info.

I found that using OpenStreetMaps (OSM), we could download the data containing all the info, so I downloaded the continent wise data from the below site.

Later I stored the information into Postgis using Osm2pgsql. I am using postgis query plugins to select the information like polygon ID and other things. But how do I get the actual polygon geometry of a point containing lat/long of type=farmland.

Example:

SELECT Polygon_geom  //random name
FROM planet_osm_polygon
WHERE
 boundary = 'farmland' AND
 AND ST_Intersects(
    ST_GeomFromText('POINT(-122.679 45.519)', 4326),
    ST_Transform(way, 4326)
)

EDIT:  I am looking for something specific like this as shown below, say I have a point on a farm and I don't know the boundary.

And using OSM (Planet-OSM) data, I get the boundary. But not what I expected (querying for landuse='farmland') doesnt return result in few cases. Does the Planet-OSM data contains the polygon info of all the farmland in specific? or it depends on the country wise data with only administrative boundaries? Is it the right way or is there any other alternate solutions?

EDIT:2
As per NMScomment got the output from overpass as follows. But my main concern still remains unanswered. How do I get an individual polygon? Now I either get a polygon with more than 3 in a single polygon or nothing (maybe OSM considers them non-farmland). Do I need to filter my query further. Please correct me if I am going in a wrong direction.


texts: Is there something like a fast access for QGIS?

texts:
When I need QGIS only to check up something, opening the QGIS application takes way to long (on my computer: 5 minutes from double cklick on the the project until ready).

I tried QGISCloud which doesn't work. Is there anything similar to just display data and maybe select different layers?

I have approximately 60 point vector layers with 30 objects each. Further a few point layers with 100-300 objects. Further, each point vector goes with a line vector (created with a query) to connect point and label. Last, there are few normal line vector layers and one google maps layer. The project itself is around 4mb.

it might work with the «Export to OpenLayers 3»-Plugin.


texts: Converting 9319106.045445956 in unknown coordinate system to WGS84?

texts:
From some API I'm receiving latlng coordinates in a stranger format: x: 9319106.045445956, y: 7020308.253131529, I don't know what it is, but my map module works only with WGS84 coordinates: lat: 53.2011871329117, lng: 83.71485471725465

What is original coords system and how to convert these coordinates to WGS84?


texts: Updating ArcGIS Online data from SQL Table

texts:
I am attempting to build an Operations Dashboard in our ArcGIS Online account using data from a SQL Server database. I’ve successfully added the data from SQL to ArcMap using the New Query tool, and published it to my AGOL account as a service. However, the data now on AGOL is static, it represents the data at the time I published it. The SQL database is constantly being updated (911 calls). Is it possible to get the data on my AGOL account to automatically update to represent the most current data stored in the SQL table?

*Note: the organization does not have an ArcGIS Server.


texts:
I was wondering what CRS to use when focusing on a large scale area for an analysis.
For example the Mediterranean Sea.

I need to find the distance of each sea port from another.
That means I need the distance from a port based, say in Greece, to another one based in Spain.

My first guess was using WGS84 ESPG:4326 and then using $length to convert degrees into meters (not the most accurate process I admit).

I can't seem to find a CRS covering the whole Med Sea.


texts: From GCS_WGS_1984 to UTM WGS 84 zone 33N

texts:
I was doing my project and I had to transform a shapefile from Bessel to gcs_wgs_1984, that went fine. Now I have to put it again to projected coordinate system.

Please could someone help me transform it from WGS_1984 to UTM wgs_84 zone 33N?


texts: Trying to set custom projection for shapefile in QGIS

texts:
I'm trying to open and use a dataset (in the form of a shapefile) (https://webarchive.iiasa.ac.at/Research/FOR/forest_cdrom/data/landscapes/landscapes.zip) with QGis. The readme of the dataset says following about projection:
Projection: Lambert Azimuthal
Units: Meters
Datum: None
Parameters:
6370997.24063 (radius of the sphere of reference)
100 0 0.000 (longitude of center of projection)
45 0 0.000 (latitude of center of projection)
0.00000 (false easting (meters))
0.00000 (false northing (meters))

I tried to create a custom CRS with the following parameters:
+proj=laea +lat_0=45 +lon_0=100 +x_0=0 +y_0=0 +ellps=GRS80 +r=6370997.24063 +units=m +no_defs
With this custom CRS applied to the dataset it looks shifted along the basemap.
Am I wrong with the projection string? (I also tried to change ellipsoid to 'sphere', without visible improvements.)


texts: Are ArcGIS Pro Pop-ups same as ArcMap HTML Popups?

texts:
I am trying to come up to speed on ArcGIS Pro Pop-ups and a piece of information that eludes me is whether reading any ArcMap documentation on HTML Popups is likely to be helpful or confusing.

Can ArcGIS Pro Pop-ups be considered to work the same way as ArcMap HTML Popups or should I forget what I have learned or read about ArcMap HTML Popups and focus solely on ArcGIS Pro Pop-ups documentation?


texts: Saving currently opened mxd and changes using ArcPy?

texts:
Here is what I have so far...

import arcpy

import arcpy.mapping as map

mxd = map.MapDocument()

mxd.save()

And yes I know that map.MapDocument requires an argument, but I want it to automatically select the currently open mxd, and overwrite it. If I can't get it to select the currently open mxd I can use:

> import arcpy
>
> import arcpy.mapping as map
>
> mxd = map.MapDocument(C:\GIS\Untitled.mxd)
>
> mxd.save()

But I would still like it to overwrite the mxd with the changes I just made to it in arcmap. Think of it like the save button except scripted.


texts: I have a web map application which displays WMS layers from GeoServer. I am looking to add a functionality which will allow users to click on multiple map objects. The application will then generate a bar chart based on the data of the selected objects and presented it to the user. Any hints or guidance on how to achieve this would be greatly appreciated.


texts: QGIS large map with high resolution details

texts:
I am trying to make a map that shows a very large region of southern Ontario but i would like to be able to zoom in on the map once it is created to see fine details such as the names of much smaller roads and streets.
When i tried to export as a pdf with a scale of 1:5000 with the appropriate map extents the resulting pdf is completely blank (although it takes awhile to save).
When i try to do a similar process in print composer the map scale seems to be different then the print composer scale.
Ive made the print composer sheet very large 20mx10m but it doesnt result in any more detail being shown, it simply is resulting in a larger map.


texts:
I want to extract the GetMap URL from WMS layers with PyQGIS. For most of the WMS
layers I can extract it from layer.source(). However this does not
work for WMS layers where GetMap and GetCapabilities is different.
Like for this layer.

The only place where I can find the GetMap URL is inside the
layer.htmlMetadata() document. From there the URL can be extracted
using regex. But every localization of QGIS (e.g. German, French ...)
has different terms for the surrounding GetMap HTML element. This
makes this approach quite tricky.

Is there any more elegant solution how to retrieve the GetMap URL from
WMS layers with PyQGIS?

PS 1: I created a workaround with OWSlib, but I am not happy with it, because the WMS needs to requested again, even though QGIS has already requested it.

PS 2: I asked the same question on the mailing list qgis-dev.


texts: Stopping arcpy for loop in cases where outputs match wildcard and get fed back to for loop list

texts:
Using a list of the feature classes matching a wildcard, I loop through 12 feature classes with a dataset. For each feature class, I make second list of all unique values within the REGION field, and then copy the feature class into mutually exclusive subsets (essentially breaking up the feature class using the REGION values). My script runs and creates all of the desired new feature classes, but I get the following error:

line 42, in  new_feature_class(),
line 28, in new_feature_class arcpy.Select_analysis(in_features, out_feature_class, where_clause),
line 90, in Select
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000670: output Output Feature Class is same as input Input Features
Failed to execute (Select).

My guess was that since the outputs are being saved into the same folder, that they are getting fed back into the for loop list, but an output can't be made because the output feature class that it would make already exists.

I'm wondering if there's a way to eliminate getting this error-- perhaps by stopping the for loop in some way?

I've copied my code below in case it's helpful.

import arcpy
from arcpy.sa import *

arcpy.env.workspace = 'C:/classifications.gdb/cy2016'
arcpy.env.overwriteOutput = True

def new_feature_class():
#loop through each feature class ending in _SIG
    wild_card = '*SIG'
    featureClassList = arcpy.ListFeatureClasses(wild_card)

    for j in featureClassList:

        # Create list out of all unique values in REGION field.
        def unique_values(table, field):
            with arcpy.da.SearchCursor(table, field) as cursor:
                return sorted({row[0] for row in cursor})

        uniqueRegionList = unique_values(j, 'REGION')

        # 3. For each feature class, create a new feature class subset for every unique region within it
        for i in uniqueRegionList:
            in_features = j
            # Make the output name: region + last 11 characters from feature class
            out_feature_class = str(i) + in_features[-11:]
            where_clause = " \"REGION\" = '%s'" % (i)

            arcpy.Select_analysis(in_features, out_feature_class, where_clause)

            # Push new subsetted feature class into next step
            in_table = out_feature_class

            # Select rows where class = 4 and clear values
            field_names = ['CLASS', 'HEIGHT']
            with arcpy.da.UpdateCursor(in_table, field_names) as cursor:
                for row in cursor:
                    if (row[0] == 4):
                        row[1] = 0
                    cursor.updateRow(row)

# Execute function 'new_feature_class'
new_feature_class()


texts: QGIS 2.14 does not work after installation

texts:
After installing for the first time qgis 2.14 on my laptop, starting qgis with the shortcuts is not working, no error messages. (Windows 7, ArcGIS Desktop already installed, 64 bit System). I have no idea, what could be wrong.
(QGIS had not been installed before on my Laptop, QGIS Standalone Installer 2.14 (64-bit))
Corinna


texts: create output line layer between two overlapping line layers

texts:
I have two street layers in QGIS. My main aim is to create an output layer of the overlapping streets (extract the overlapping portions from the red streets via blue intersect layer). The file size of the red streets is too large, therefore I want to select just the streets, which are within the map extent for my intersection.  The output should also contain the attribut information of the red streets.

I copied a code in which I could intersect both layers, but I could not fix the selection within the map extent, and I could not create an output layer with the attribut information. Below you can see my code and an image:

registry = QgsMapLayerRegistry.instance();
blue_street = "blue_streets";
red_street = "red_streets";

route_blue = registry.mapLayersByName( blue_street ) [0];
streets_red = registry.mapLayersByName( red_street ) [0];
extent = iface.mapCanvas().extent()

selections= []
for f in route_blue.getFeatures():
    for a in streets_red.getFeatures():
        if a.geometry().intersects(f.geometry()):
            intersection = a.geometry().intersection(f.geometry())
            print intersection.exportToWkt()
            selections.append( f.id() )
        break


texts: In GMT after adding scale bar , psxy plots are disappears

texts:
Here My script is

!/bin/bash

LATLON=30/130/-40/20

PROJ=M6i

STASYM=t0.1

STACOLOR=255/0/0

FNAME=map.ps

EPISYM=c0.1

EPICOLOR=255/255/0

gmt grdimage ../../ETOPO1.grd -R${LATLON} -J${PROJ} -C../../sediment.cpt -P -B30 -V -K -Y10 > ${FNAME}  << EOF

EOF

gmt psscale -C../../sediment.cpt -D7.7/-1.0/15/0.4h -P -O -B2000:"Bathymetry & Topography relief":/:meter: -I -V   >> ${FNAME}  << EOF

EOF

#MAP GREAT CIRCLE RAY PATHS USING PROJECT
gmt psxy ./RAYPATH -P -J${PROJ} -R${LATLON} -W0.1 -O -: -V -K >> ${FNAME} << EOF

EOF

#PLOT EPICENTER LOCATIONS

gmt psxy ./EPICENTRE -P -J${PROJ} -R${LATLON} -O -: -S${EPISYM} -W0.2 -G${EPICOLOR} -V -K >> ${FNAME} << EOF

EOF

#PLOT STATION LOCATIONS

gmt psxy ./STATION -P -J${PROJ} -R${LATLON} -O -: -S${STASYM} -W0.2 -G${STACOLOR} -V   >> ${FNAME} << EOF

EOF


texts: Earth Engine: Exported image showing up as white image

texts:
For Landsat, I am trying to export a rectangular geometry that I drew for the California Kingfire (prefire) as a GeoTIFF. I am able to successfully export it to my drive, but the image is white. It seems that I might be only looking at one band (maybe?). How can I show the full colored image?
This is my code in the Earth Engine code editor:
var dataset = ee.ImageCollection('USDA/NAIP/DOQQ')
                  .filter(ee.Filter.date('2012-01-01', '2013-12-31'));
var trueColor = dataset.select(['R', 'G', 'B']);
var trueColorVis = {
  min: 0.0,
  max: 255.0,
};
Map.addLayer(trueColor, trueColorVis, 'True Color');

// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: dataset.mean(),
  description: 'Kingfire_RGB_pre',
  scale: 30,
  region: geometry,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

Any ideas?


texts: QGIS - WHERE SQL within a processing.runalg statement for v.surf.rst in python console

texts:
I'm creating lots of raster files with grass:v.surf.rst so rather than using the graphical interface I wanted to use the python console (and potentially put in in a script at some point).  However when I try and run a script copied from the history I get a syntax error if I include the WHERE SQL syntax.  If I remove "x" <> 3.54 and just leave empty double quotes it works fine.
An example of the code:

processing.runalg("grass:v.surf.rst","file://drive/filename.csv?type=csv&xField=x&yField=y&spatialIndex=no&subsetIndex=no&watchFile=no",""x" <> 3.54",None,"DD275",40,40,300,0.001,2.5,1,0,0,False,False,"1.5,5.25,20.5,32.5",0.1,-1,0.0001,"//drive/filename.tif",None,None,None,None,None)

Is there a different process for the where statement in the console?


texts: Incorporating SQL Statements in ModelBuilder?

texts:
I have a point data set where each point has a unique operator associated with it.  I'd like to select by attribute = operator name and generate a new shapefile for each operator containing only the points associated with that operator name.  Doing this by hand for each operator is simple but I've been having trouble trying to build a model in ModelBuilder that will automate it.

I've tried using the Select Layer by Attribute tool but I can't figure out how to iterate it such that the SQL expression will run through each operator name.  ModelBuilder will not let me attach any inputs onto the expression term.  The gist of what I think I'm trying to do here:


texts:
Does anyone knows if there exists a caching system for views or queries in postgis?
I have a system that gives information through a webpage using a postgis view as backend. The problem is that, sometimes the view is big enough that takes 5 seconds per user to serve the info. I know there's this system named geocache that works very good for raster data using wms protocol. I was wondering if there's a similar tool for vector / table data.


texts: Open source absolute atmospheric corrections for Sentinel-2

texts:
I am searching for open source absolute atmospheric correction processors for Sentinel-2. There are many of proprietary ones (FLAASH, QUACS, MACCS) which I cannot use for my study.

Beside DOS1 (QGIS SCP plugin), Sen2Cor (SNAP) and the upcoming new 6S/ATCORR (GRASS), are there any other absolute atmospheric correction processors for Sentinel-2 imagery available in open source software packages?


texts: I have a CSV file with geocoded data that I would like to display in QGIS.  How do I open this file?  When I go to add CSV layer it asks for X and Y co-ordinates which I do not have as the data is geocoded?


texts: Get maximum pixel value within numeric range

texts:
I am working with MODIS Global Primary Production (GPP) rasters and I am trying to pick the max GPP pixel value in a series of 15 dates for a study area and create a new raster layer. MODIS GPP rasters contain values greater than 30,000 that are either NoData, unclassified or other categories that do not correspond to GPP values. I am actually looking for the max value between 0 and 30,000.

gdal_calc.py will only take 2 rasters as input and while SAGA accepts more than two rasters as input, it does not let specify a range for values to be taken into account.

Any suggestions how this can be done?


texts:
I am trying to make my own vector tiles and I am using Mapbox Studio Classic for that purpose.
When I export my dataset, Mapbox Studio will try to figure out the extend of my data and those queries take very long.

The query it's performing is like this:

SELECT
  ST_XMin(ext),
  ST_YMin(ext),
  ST_XMax(ext),
  ST_YMax(ext)
FROM (SELECT ST_Extent(geom) as ext from cont) as tmp`

And my cont table is quite huge with a lot of short (non of them is long) contour lines in it.
It has an index on the geom column (GIST).

The query is taking a LOT of time (several minutes) on our beefed up machine.

Is there a way to improve the performance of that query?

EDIT:

EXPLAIN (ANALYZE, BUFFERS) generates:

Subquery Scan on tmp  (cost=13016236.00..13016236.04 rows=1 width=65) (actual time=2245959.393..2245959.396 rows=1 loops=1)
  Buffers: shared hit=39596394 read=21929928
  ->  Aggregate  (cost=13016236.00..13016236.01 rows=1 width=5190) (actual time=2245952.748..2245952.748 rows=1 loops=1)
        Buffers: shared hit=39596394 read=21929928
        ->  Seq Scan on cont  (cost=0.00..12830974.00 rows=74104800 width=5190) (actual time=0.014..104310.564 rows=74119236 loops=1)
              Buffers: shared hit=3343020 read=8746906
Planning time: 0.086 ms
Execution time: 2245959.495 ms


texts:
When running my Python program, the following error came out:

ERROR 1: NUMPY driver was compiled against GDAL 2.1, but the current library version is 2.2

What can I do to fix this?


texts:
I have a dataset that I have pulled from a MaxSea 3d database.

I am trying to figure out what projection is used, so I can transform and use the data. I have tried following this guide http://support.esri.com/en/knowledgebase/techarticles/detail/24893, for a small sample of data (300k points). Firstly, the fields are named lat and lon, indicating a Geographic Coordinate System (datum), but the data given looks more like it is in a Projected Coordinate System, becuase the numbers are in 8-9 digits range.

The data on file is collected around the Faroe Islands, and if the data is in Geographic Coordinates, it should be in or around 62°N and 7°W ± a few degrees, so I can conclude that the data is not in any format where the leading digits are degrees. (not decimal degrees or DDMMSS.ss, as suggested here).

That leaves me with either a proprietary, unknown projection (possible!), or it is in a Projected Coordinate system. I keep on following the Esri guide, to figure out what Projection is used, and this is where I need help, because the Esri guides only covers the U.S.

  If the data is in the United States and shows an extent in which the coordinates to the left of the decimal are 6, 7 or 8 digits, the data is probably projected to the State Plane or UTM coordinate systems.

Here I see that there are 9 digits (!) left to the decimal point, and when I try to add some reference data, I find that my data is outside the extent of the earth.

What would you try, I you were asked to figure out the projection of the data?

Here is a sample of the data (50k points). The fields are: id, level, lat, lon, weightdepth_cmap, valuedepth_cmap, weightdepth, valuedepth. The fields of interest are lat, lon.


texts: QGIS hyperlink to images using EVIS and relative path

texts:
Just wondering if anyone has experience creating relative paths in the EVIS plugin for JPEGs.
Basically got a set of points that id like to have linking to some JPEGs. They are all stored
in the same database but not in the same folder.
It's been a few years since I've had to use relative paths and cant quite seem to track down the reason why its not working. Is there a general rule for the end portion in the image path column? e.g., just the folder the images are in or it and a set number of layers above?


texts: PostGIS ST_AsTIFF function raises memory exception?

texts:
We have a satellite image loaded into PostGIS.  We've calculated a reflectance raster from that image and stored the results as 32BF pixels.  There are 54 1000x1000 tiles in the reflectance raster.

I want to create a GeoTIFF of that reflectance raster.  Individual tiles work fine.  But using an ST_Union on the reflectance raster column raises an exception when executing the query below from Python:

SELECT ST_AsTIFF(ST_Union(refl_rast), 'GTiff') gtiff
FROM sat_rast_calc

InternalError: invalid memory alloc request size 1344128795

I know I can export the tiles individually and merge them with gdal_merge.py afterwards.  I'd prefer to export this GTiff in one step though.

Would any of these be likely to help?

Using a different pixel type
Passing in another ST_AsGTIFF option
Increasing a PostGRES memory parameter
Increasing a Linux kernel parameter
Doing something else


texts: ZSTD support for GDAL GeoTIFFs with UbuntuGIS

texts:
I have some GeoTIFFs that are compressed with the ZSTD standard. I'm trying to serve those on a VM that has UbuntuGIS installed via the UbuntuGIS PPA: https://launchpad.net/~ubuntugis/+archive/ubuntugis-unstable on a base Ubuntu 18.04 image. The GDAL version on that PPA is 3.0.4.

When I try to gdal.OpenEx this file on the VM I get a ERROR 1: Cannot open TIFF file due to missing codec. error because I am missing the ZSTD codec I assume.

Is there a way to add ZSTD support for GDAL in UbuntuGIS? Perhaps through an apt install that I'm missing?


texts: Python To PostgreSQL - ST_GeometryFromText and NaN

texts:
I am trying to import a shapefile into Postgre using python and psycopg2. I have run into a few issues.

I am using ST_GeometryFromText() of a WKT string.

First, the shape@WKT returns MULTIPOLYGON ZM at the beginning of the string. I believe that it is only supposed to be MULTIPOLYGON correct?

I fixed this by doing a replace on the string but this does not seem like what should be done. What do I do to make this correct?

Second, the shape@WKT contains X, Y, 0, NAN as values. I am expecting only 2 values, x and y. Why am I getting 0 and NAN?

with arcpy.da.SearchCursor(input_fc, ["FID", "Shape@WKT"]) as da_cursor:
    for row in da_cursor:
        geoms = row[1].replace(' ZM ','')
        print geoms
        pg_cursor.execute("INSERT INTO %s (fid_p, geom) VALUES (%s, ST_GeometryFromText('%s', %s))" %(dataset_name, int(row[0]), geoms, int(coordinate_system)))
pg_connection.commit()

RESULT of print geoms:

MULTIPOLYGON ZM (((440185.43178632273 6350644.499023444 0 NAN, 440189.09859540447 6350575.940429694 0 NAN, .........

PYTHON ERROR:

InternalError: parse error - invalid geometry
HINT:  "...85.43178632273 6350644.499023444 0 NA" <-- parse error at position 60 within geometry


texts: Somehow I've created two projects with duplicate data in QGIS. How can I delete one of them?


texts: Beautiful MapServer maps and their map files

texts:
Looking for

beautiful MapServer produced maps

and

their associated map files to get insights from it.

beautiful is subjective,

quick definition: pleasant to look at such as maps from http://mapbox.com/maps/

Any ideas where to get them ?


texts:
I am trying to import OSM data into Spatialite.

The data is a query from overpass called features.osm

I'm using spatialite 3.0, and following examples, I have created a new database with:

$ spatialite features.sqlite

then I import the data using

$ spatialite_osm_net -o features.osm -d features.sqlite -T features

This results in the following output:

SQLite version: 3.7.9
SpatiaLite version: 3.0.0-beta

Loading OSM nodes ... wait please ...
    Loaded 11 OSM nodes

Verifying OSM ways ... wait please ...
    Verified 0 OSM ways

Disambiguating OSM nodes ... wait please ...
    No duplicate OSM nodes found - fine ...

Loading network ARCs ... wait please ...
    Loaded 0 network ARCs

Dropping temporary table 'osm_tmp_nodes' ... wait please ...
    Dropped table 'osm_tmp_nodes'

Dropping index 'from_to' ... wait please ...
    Dropped index 'from_to'

VACUUMing the DB ... wait please ...
    All done: OSM graph was successfully loaded

However when I open the database and view the contents of 'features' table, its empty. What am I missing here? Perhaps I'm trying to fit a square peg in a round hole?

What I want, is to be able to store osm data in the sqlite database and execute queries to retrieve nodes/ways/relations for a specific bounds (just like my original url request).


texts: Create a polygon from an array of coordinates in arcgis javascript api

texts: Well I am building a GIS web app using ArcGIS Javascript API , I want to create a user interface in which the user provides an array of coordinates and I have to get those coordinates and create a polygon on the map from this array of coordinates. Do you have any idea?


texts: Populating field with the name of the input file

texts:
I'm using modelbuilder to import a dbf (selected by the user) to a geodatabase and turn it into a feature class for further processing.

I would like to be able to populate a field in the imported table with the original file name of the source data.

This is section of my model where it imports the dbf (named Input Rows by the model) and adds the field but I can't find a way of capturing the input file name to use in Field Calculator

Edit: I've looked at this, but it seems to be talking about adding a field based on the feature class/shape file name. My table will have a generic name, I want to be able to use the name of the file that is being selected to populate it.

The geodatabase is a template that will be re-used by lots of users in different locations with differently named dbfs.


texts:
I am putting together a compilation of 1/3 arc second, and 1/9 arc second DEM data for essentially all of southeastern Oregon in ArcGIS. This surface would be used for geomorphological and geological modelling of various types, and to show areas and points where previous archaeological and geological research has been done.

The question is: what are the best practices one would use to put together such a large mosaic of high resolution DEM data? Is there any recommended literature I should check out? Optimally, this would all be in one map, but I worry about it being unwieldy and taxing on memory, especially when analyses are being run. Perhaps it should be broken down into subregions of some type?


texts: Batch Clipping Single Raster with Multiple Shapefiles in ArcPy?

texts:
I have a single raster file and about 40 small polygons shapes. I would like to clip the raster by the extent of all these shapes, and create 40 resulting .img files.

I know I can just do this by right-clicking the clip tool and selecting batch process.... but I would rather have a script because I will need to do this multiple times over the next few months.

I am inexperienced in python and would like some help with my script. I am not using geodatabases, just a simple input and output folder. I was looking at a similar post which has the below code... but the rectangle number option in arcpy.Clip_management requires me to manually enter every single rectangle number which defeats the purpose for me.

How can I modify this to manually clip a single raster using all 40 polygon files?

I am using ArcGIS 10.2.2.

import arcpy

rasterlist = arcpy.ListRasters("your_workspace")
for raster in rasterlist:
for i in range(40):
    arcpy.MakeFeatureLayer_management("your_shapefile", "layer" + str(i), ' "FID" = ' + str(i)) #create a layer with only polygon i
    arcpy.Clip_management(raster, "#", raster[:-4] + "clip" + str(i) +" .img","layer" + str(i), "0", "ClippingGeometry") #clip based on layer, clipping geometry will use the polygon extent only


texts: Leaving Image Objects that are Uncertain to be a Particular Class as 'Unclassified' in eCognition

texts: In classification of the imagery, is there a way of leaving certain image objects unclassified? Say for example, if the classifier is uncertain what class to classify the image object as, does eCognition have functionality to leave the image objects unclassified? I'm not sure if that is clear, but let me know if you all have any ideas. I've looked into hierarchical classification, however, I am using a ML classifier, RF, to do the classification, which is separate from the hierarchical classification algorithm. The eCognition Community site has been down for several weeks now, and there is sparse information across the web, as well as the User Guide, which is not very helpful. Thanks in advance.


texts: geoxml3 parseKMLString function not working?

texts:
I am unable to use Parser.parseKMLString().
Its giving me a error saying Parser.parseKMLString() is not a function.
What i am trying to do is, read the KML file on the server and return the entire file as a json string which is then passed to the parseKMLString function.
But i am getting the above error.


texts: Searching for NULL values in feature class table

texts:
I am looking for a way which I can search whole my feature class in ArcMap to find if has any NULL values on it.
I know how can I do this with "Select By Attribute", but I am looking for a faster way to do this, like with a Python script.


texts: Extracted by mask raster file with evapotranspiration data not displayed

texts: Should extract a part of raster file containing evapotranspiration data with a polygon.The polygon is that of a barren land which was made with tool ,extract by attributes from my Land use classified raster file and then converted to polygon with Raster to polygon conversion tool.But after doing extract by mask with the same polygon, the output raster is not displayed.I am using Arcmap 10.4 software.


texts: How to output the number of channels in a .JP2

texts:
I am looking for any tool (or python code snippet) that can loop through a file folder of .jp2 images and determine which have 3 bands and which have 4 bands.

I am having trouble creating a raster dataset as some of them have 3 bands and others have 4. I need to isolate the band count and replace the images with only 3 bands. I am looking to output a text file that has the name of the file which has 3 bands.

Any thoughts?


texts: Creating TIN from UAV images or Google Earth images?

texts:
I'm starting using ArcGis and I would like to create a TIN from UAV images.

What are the steps that I have to follow?
I think:

Mosaic the images, merge o join. (I have 40 aereal pictures in .JPG)
Georeference the big joined picture with GPS points (I may have these points, but in other cases not).
Do something else.
Create the TIN or xyz cloud points to export to AutoCAD.

How can I do this?
Do I need other tools?


texts: Is there a way to "wrap on character" in the identify result box?

texts:
I have lists of plant species recorded at different sites (polygons). I have listed the species with their coverage in a seperate field. species seperated by semicolon. When showing the information of the polygon "identify results" the list is cut (see picture).

Is there a way to wrap the list on character to get a list in vertical column alignment. Best would be after every second semicolon.


texts:
í´ve moved several Points from my Survey file. Now i want to recalculate the new coordinates of those moved Points. And then Export it to Excel.
Does anybody have a quick answer?


texts: Convert all .lyr in an .mxd into kmz using arcpy

texts:
I wrote a simple script supposed to convert all layers in a mxd into kmz and it doesn't work.
Kmz are being created but they're not what they should be.

import arcpy

arcpy.env.workspace = pathworkspace
mxd = arcpy.mapping.MapDocument(pathmxd)

allLayers = arcpy.mapping.ListLayers(mxd)
for layer in allLayers:
    outputkmz = layer.name[:-4]+".kmz"
    arcpy.LayerToKML_conversion(layer,outputkmz)

If i try to open myself in ArcMap the kmz created by the script, I get this error message.

Some points are being displayed but in a strange way, like if they were aligned while it's not the case.
If I open them with Google Earth, it just doesn't display anything.

Other information:
- I'm using I'm using ArcMap 10.5.

I'm using PythonWin
If I try the script with a .mxd where there's only one layer in the ToC, it's working. As soon as I add another layer, it doesn't work anymore.
If I try the conversion manually with the toolbox in ArcMap (convert layer to kml), it's working, so I guess it means the problem doesn't come from my .lyr files.
Doesn't seem to be a problem of coordinates system

This question is related to my other (more global) question here: Convert Shp to KMZ displaying labels using arcpy?


texts: I have a map georeferenced in georeferecer plugin in QGIS. I have loaded it in qgis. It shows strangely (wrong aspect ratio than I expect). I am changing projection (bottom right "OTF" projection) from "Gauss Kruger zone4, to WGS).  This first results to change of position of the window to some ridiculous coordinates, after I find my original coordinates I see matched points but I see no raster. Why?


texts: Reprojecting to same srs shifting the layer

texts:
I am using ogr2ogr to reproject a shapefile. Here is what I have done.

I have a shapefile A.shp. I used the following command to create B.shp:
ogr2ogr -t_srs EPSG:3857 B.shp A.shp
Now I reprojected B.shp again and created C.shp
ogr2ogr -t_srs EPSG:3857 C.shp B.shp

Since B and C have same data and they are in same projection, I think they should completely overlap. However the do not. In the image bellow the Red one is A.shp and the other one is B.shp.

Is there a good reason why this can happen? This is happening to all the shapefiles not just this one.


texts: QGIS: Update geometry via geom_from_wkt('POINT(X,Y)') - how to correctly parse string columns?

texts:
The columns containing X and Y strings look like

-1234,56

So far i tried:

geom_from_wkt(concat('POINT(',replace("wbr_komplett.X",',','.'),',', replace("wbr_komplett.Y" ,',','.'),' )'))

and also:

geom_from_wkt('POINT(
'|| replace("wbr_komplett.X",',','.') ||','||replace("wbr_komplett.Y" ,',','.') || ' )')

Both to no avail. After saving the edits and closing the session the desired features still don't have a valid geometry. The layer in question is a Geopackage. As far as i can tell both queries are correct - also they wouldn't be executable if they were not.


texts: How to fix encoding problem in shapefile importing?

texts:
I tried to import a shapefile to postgresql. I used PostGIS shapefile import/export manager. Earlier I have done same thing but did not face any problem. But this time I am facing problem with encoding. First I tried to import into an exsisting database with UTF-8 encoding. It showed following error:

Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Unable to convert data value to UTF-8 (iconv reports "Invalid argument").
Current encoding is "UTF-8". Try "LATIN1" (Western European), or one of the
values described at http://www.postgresql.org/docs/current/static/multibyte.html.

Shapefile import failed.

So, I created another database with 'LATIN1' encoding with collation C and character type C. Then tried to set client encoding to 'LATIN1' by executing the following query.

set client_encoding='LATIN1';

Now tried to import. But still it is showing client encoding is UTF-8 and import fails. Screen shot is attached. What to do?


texts: Refreshing ArcIMS service with ArcObjects

texts:
Does anyone know if it is possible to refresh an ArcIMS service with ArcObjects and if so do you have a nice code snippet? :)

I think it might have something to do with interface IIMSAxlRequest - but not sure how to use it...


texts: Use single sign-on username in Web AppBuilder filter widget?

texts:
I have an Portal web map:

Details:

The map has a Web AppBuilder filter widget
It uses single sign-on (SSO) authentication

The map has a work order layer that is integrated from a Workorder Management System (WMS). Each work order record in the layer has a work order username value (example: JDOE).

I would like to configure the map so that it:

Grabs the SSO username that it used to authenticate the web map
Filters the work orders based on the username

Is it possible do this--to filter the layer based on the SSO username?


texts: Failed to connect to database. Not a valid file name (ArcMap, personal geodatabase)

texts:
Windows 10, ArcMap 10.8, local parent directory, personal geodatabase
I tried solutions suggested for this problem but they don't work. I have several similar ArcMap projects for different regions. In some projects layers are not accessible and when I try to change the data source an error pops up: "Failed to connect to database. An underlying database error occurred. Not a valid file name".
I tried changing the database file name, tried to change the directory of the file, tried to run AccessDatabaseEngine but nothing works. Some solutions suggested deleting the file lock but I don't have it in the file directory. When I open the DB in Access I can see the attributes, it doesn't suggest me to fix the DB (I don't have licensed Office but other ArcMap projects still work).
Those projects that DO work are dated this year and those that DON'T are dated last year so I thought there could be a compatibility problem. Unfortunately, because I have not created those projects I don't know which ArcGIS version was used.
How can I fix this issue?


texts:
I would like to use color scale of my raster data in my leaflet web application exactly same way as it is in this video https://youtu.be/l47ABvDLiq8?t=51s .

Does somebody know which service is used for raster data publishing? Is it basic geoserver? Or is there used some another wms service?


texts: hdf to gtiff, reprojecting, cropping with a shapefile in batch using R

texts:
I have some raster data from MODIS, with hdf format.
I want to convert them into GeoTIFF, reproject them and crop them based on my shapefile.
I was able to do it in GDAL, in batch mode, with the following code, but I don't know how to do it in R.

gdal_translate -of GTiff input.hdf output.tif
gdalwarp -t_srs "+proj=utm +zone=16 +datum=WGS84" input.tif output.tif
gdalwarp -cutline tl_2013_18_cousub.shp -crop_to_cutline -dstalpha input.tif OUTPUT.tif

I use the following code in R, but it's just for the converting to GeoTIFF and also it does not provide the result with TIFF format. (I don't know why)

    library(gdalUtils)

files <- dir(pattern = ".hdf")
for (i in 1:length(files)){
  sds <- get_subdatasets(files[i]) # the input files have different subdatatsets, and I need only the first one
  gdal_translate(sds[1], dst_dataset = paste0("new_",basename(files[i])), of = "GTiff")
}

based on your comment, I wrote the following code though it has errors,
can anyone help me with the issues in comments?

   library(gdalUtils)
        files <- list.files(path=".../test", pattern="*.hdf", full.names=T, recursive=FALSE)
    for (i in 1:length(files)){
    sds <- get_subdatasets(files[i])
    gdal_translate(sds[1], dst_dataset = paste0(basename(files[i]),".tif"), of = "GTiff")
    }
#first issue: the file name would be like "filename.hdf.tif" how can I remove "hdf" part to make it like "filename.tif"

    shp <- shapefile("tl_2013_18_cousub.shp")

#second issue: I had to write another loop, can it be merged with previous loop to save time in running process?
    files <- list.files(path=".....", pattern="*.tif", full.names=T, recursive=FALSE)
    for (i in 1:length(files)){
    sds <- raster(files[i])
    xp <- projectRaster(sds, crs = "+proj=utm +zone=16 +datum=WGS84")
    writeRaster(xp, file = paste0("prj_",basename(files[i])))
# third issue: the saved file has different resolution! how can I maintain the original resolution

    cr <- crop(xp, shp)
    writeRaster(cr, file = paste0("crp_",basename(files[i])))
# Forth issue: it guves me error: "Error in .local(x, y, ...) : extents do not overlap"
    }


texts: Joining\relating\connecting 4 feature class together?

texts:
How can I join\relate\connect  communities\UMP\BuiltUP\Admin boundaries?
Could you please consider the following scenario:

I have a feature class with a name “community “ which I consider it the main layer in may work, so I relate it to more than one table and it has many information about the community.

Then I have three other feature classes each one is stored in a deferent file geodatabase ( built-up area, UMP urban master plan, and Administrative boundary).

Now how could I join\relate\connect the three feature classes (built-up area, UMP, and Administrative boundary) to the (community) feature class so that when I identify the community point I can get the whole results from all feature class?


texts:
Does Grass GIS consider a 3D set as a single XY location, with a Z (elevation), and a property P.  I have a dataset with vertical data points (stacked) in XY coordinates.

Each XY location has a Z and a Property. therefore i have as many XY as P.

I am wondering if I should slice the dataset into individual elevations (one data set for each Z), so that I can import and interpolate.


texts:
we have built a model with pythn/GRASS that calls the r.out.mat module. The area is quite large 25,000x25,000 cells. The machine we are running it on is Windows 7 machine with 16GB of RAM.

After a while the script aborts complaining about an out of memory issues. Now, I would not have expected to have such issues with so much RAM even though the manual page for r.out.mat does say that the module loads the whole map into RAM so it requires a huge amount of RAM. Any idea on how we can optimise the script avoid such problems? If this is a limitation of the module, I was wondering whether r.out.mat could be used to read the map area in patches, dump out the results in single files and then some other module patches all the r.out.mat results together. Any suggestions on how we can overcome this issue?

Thanks
Antonio


texts:
When I work in QGIS and try to fill gaps between polygons creating another one over them (project snapping settings are selected as Avoid intersection), I hope "fill" this gap with just one polygon with the same shape, but sometimes when I create this polygon I create a multipart polygon, which contains the principal polygon and some segments of 2-3 vertex, most of them, between other polygons created and where there are not gaps.
Why does it happen, if snapping setting must avoid it?


texts:
I have been looking at D3 and following along with their basic tutorials. In the documentation you can use json files to create graphs, however I am confused on how to make the connection from my cartodb query,  to the input of D3. To clarify I have made a successful query to my cartodb table and returned those results from there I need make a connection from my returned json query to D3.

  var sql = new cartodb.SQL({ user: 'ceigis' });
sql.execute("SELECT status, COUNT(status) FROM deldot_rail_crossing GROUP BY status")
  .done(function(data) {
    for (var i = 0; i < data.total_rows; i++) {
      $('#content').append("
" + data.rows[i].status + "" +  data.rows[i].count+ "

"); ///This was just to test if my json query was correct
     //Making the Connection from my query to D3 is where I am struggling

        d3.json("data.json", function(json) {

I am ultimately trying to create a simple pie graph of how many records fall within 4 different statuses.

UPDATE: I was able to get a working pie chart. One of my problems was passing json object instead of an array.code posted below.. hope this helps out others

function init(){
var w = 1000;
var h = 1000;
var r = 200;
var color = d3.scale.ordinal()
    .range(["#98abc5", "#8a89a6", "#7b6888", "#6b486b", "#a05d56", "#d0743c"]);

var sql = new cartodb.SQL({ user: 'igis' });
sql.execute("SELECT status, COUNT(status) FROM rail_ing GROUP BY status")
        .done(function(data) {
            for (var i = 0; i < data.total_rows; i++) {
                var data = data.rows;
            }
        var canvas = d3.select("body").append("svg")
            .attr("width", w)
            .attr("height", h);

        var group= canvas.append("g")
            .attr("transform", "translate(300,300)");

        var arc = d3.svg.arc()
            .outerRadius(r)
            .innerRadius(0);

        var pie = d3.layout.pie()
            .value(function(d) { return d.count; });

        var arcs = group.selectAll(".arc")
            .data(pie(data))
            .enter()
            .append("g")
            .attr("class", "arc");

        arcs.append("path")
        .attr("d",arc)
        .attr("fill", function (d) { return color (d.data.status); });

        arcs.append("text")
        .attr("transform", function (d) { return "translate(" +
                    arc.centroid(d) + ")";})
            .attr("text-anchor", "middle")
            .attr("font-size", "1em")
            .text (function (d) {return d.data.status;});

        });
}


texts: Crop a raster layer given lat/lon ranges in QGIS 3.2?

texts: I uploaded a landsat raster layer in QGIS 3.2 and I need to crop it. I noticed that for every pixel QGIS knows the lat/lon value (as expected). Since each pixel has its associated lat/lon coordinate, is there a way to specify a lat and lon range and crop the raster layer based on that range?


texts: Trying to export shapefile attributetable to JSON using GeoPandas without the geometry

texts:
I am using GeoPandas' read_file function to read a shapefile
and to_file function to save to GeoJSON.
This is working fast and good.
But I need to save only the attribute table without the geometry
or with empty geometry.
How can I do that?
import geopandas as gpd
import pandas as pd
import fiona

def main():
filepath = "D:\\DATA\layers1\\TTL_TRANSPORT.shp"
destination = "D:\\DATA\\jsonfiles\sample.json"

gdf = gpd.read_file(filepath)

df = pd.DataFrame(gdf)
columns_without_geom = gdf.columns

print( columns_without_geom)

if 'geometry' in columns_without_geom:
   columns_without_geom.drop('geometry')

gdf = gdf[columns_without_geom]
df = pd.DataFrame(gdf)
df.to_json(destination)
df.to_json(destination)


texts:
I am trying to listen to right mouse click , i tried "rightclick" as in documentation says , but it doesn't work.

Can any one say why it doesn't work ???

map.on('rightclick', function (evt) {
            var feature = map.forEachFeatureAtPixel(evt.pixel, function (feature) {
                if (typeof (feature.getId()) != "undefined") {
                    return feature;
                }
            });

            if (feature) {
                var registeredObjectId = -1;

                if (typeof MapObjects !== "undefined") {
                    for (var i = 0; i < MapObjects.ObjectResult.length; i++) {
                        if (feature.p.name === MapObjects.ObjectResult[i].Id) {
                            registeredObjectId = i;
                        }
                    }
                }

                var geometry = feature.getGeometry();
                var coord = geometry.getCoordinates();

                if (registeredObjectId === -1) {
                    popup.setPosition(evt.coordinate);
                    $(element).popover({
                        'placement': 'top',
                        'html': true,
                        'content': "Delete"
                    });
                    $(element).popover('show');
                    $(element2).popover('destroy');
                }
                else {
                    popup2.setPosition(evt.coordinate);
                    console.log(MapObjects.ObjectResult[registeredObjectId].Name);
                    console.log(MapObjects.ObjectResult[registeredObjectId].Value);
                    $(element2).popover({
                        'placement': 'top',
                        'html': true,
                        'content': "" + MapObjects.ObjectResult[registeredObjectId].Name + "
" + MapObjects.ObjectResult[registeredObjectId].Value + "
Delete"
                    });
                    $(element2).popover('show');
                    $(element).popover('destroy');
                }

            } else {
                $(element).popover('destroy');
                $(element2).popover('destroy');
            }
        });


texts: PgRouting dijkstra returning null cost values

texts:
I have been following the tutorial at https://anitagraser.com/2011/02/12/drive-time-isochrones/ in an attempt to create an isochrone map showing drive times to the nearest medical practice in my region (modifying my process where necessary for PgRouting v2.5).

Whilst calculating the drive time needed to reach each of my source nodes (in this case medical practices) from all other nodes, I found that all the cost values were returning as null except for the row following the source node (in the case of the screenshot, I would get a non-null result for row 3229 only).

Is there an obvious error in my code that would result in these null values?


texts: How does "Minimum Size of Exterior Watershed Basin" work in r.watershed command in qgis 2.18.12?

texts: I have delineated the watersheds from DEM using r.watershed command. I have specified the Minimum Size of Exterior Watershed Basin as 1500. I have read the description of r.watershed as well. I am not able to understand why there is so much variation in size of watersheds. There are tiny watersheds even of the size of one cell. How can such small watersheds form after specifying criteria as big as 1500?


texts: Large point data,not visible/queryable on zoom in

texts:
Using ArcMap. I have around 15 million points layer that I can only see them when I click on 'Zoom to layer'.
When I zoom in the points disappear, I cannot even get info of these (existing) points drawing a rectangle ('select by polygon' gives 0 records).
I have no min/max 'visible scale range' setup.
I am wondering if is a performance problem, but the software does not get frozen on zooming in, and the drawing icon looks fine, but non-active.


texts: Making one-to-many polygon interactive map using ArcGIS Pro

texts:
I have a polygon shapefile of all war zone. Each war zone has a list of allies and rivals. In which case it is a one-to-many relationship.
The goal is to make an interactive map that will allow you to click on a polygon and it will highlight the allies and rivals in two different colours.
Is it possible to do this in ArcGIS Pro without coding?


texts: How to make ArcGIS Online contents sharing level funnel through Organization Admins?

texts:
Is there any way within ArcGIS Online to notify admins of items as they are shared publicly within an organization?
We are trying to moderate our content more to make sure our users don't share things publicly that they shouldn't.
Is there a way to put such items into a queue so that admins can approve and then share publicly?
I know we can search for all items shared publicly, but it would be nice to have tools to make an admin-approval workflow for publicly shared content.
Any ideas or leads to custom scripts/tools, etc?


texts: Using SqlGeometry.MakeValid to get a counter-clockwise polygon in SQL Server

texts:
The data I've loaded from a shapefile into a geometry column in SQL Server (using NetTopologySuite) is not producing valid GeoJSON when I export it because it doesn't follow the right-hand rule.

I thought this would be easy to fix with the MakeValid() function, but it's not correcting the issue.

I tried with a simplified geometry as below:

declare @geom geometry
select @geom = geometry::STGeomFromText('POLYGON ((
0 0,
0 1,
2 1,
2 0,
0 0
))', 4326)

select @geom, @geom.IsValidDetailed()

SQL Server says it's valid (24400: Valid), so presumably that's why it's not fixing it in the MakeValid function.

What's the best way to resolve this? I've seen an option to use:
UPDATE table SET geom = geom.STUnion(geom.STStartPoint()) which works for simple polygons, but I'm not sure how well it will work for more complex geometries.


texts: Saving a calibrated route in QGIS for "future" use

texts: Can a Calibrated Vector (after using LRS) be saved with M values, so that the file does not have to be calibrated every time is used? We're currently working on State road data, and it has many records both on roads, a swell as the measure points. Is there a way that once the vectores have been calibrated, these can be saved so that we dont have to go through the process of calibrating every time we use it?


texts: Where can I find Hawaii public transit data?

texts: I am looking for public transportation routes and stops for Hawaii public transit.  So far I have found a shapefile for TheBus routes and stops for the island of Oahu from the state of Hawaii Office of Planning, but it did not include bus/rail files for the other islands.


texts: Quickly plot 600+ Post Codes in ArcView 9.3

texts:
I have 613 post codes covering the UK, stored as an Excel file. I need to plot them - preferably all at once - in ArcMap (ArcView, 9.3). The post codes are 5 figures (e.g. BT19 6).

Is there a quick way of doing this? Currently I'm using Arc's search tool to find the post code, plot a point, and convert those points to a shapefile. However, with over 600 to plot, this will take an age.

I don't really have any experience of Geocoding, but I'm assuming this is the way to go?

Any advice is greatly appreciated!


texts: Replace empty values with the next available one using FME

texts:
how can be replaced empty attributes with the next available values ?

I tried

  NullAttributeMapper

but somehow I dont get desired results!

What I have:

attribute1, , ,attribute4
 ,attribute2,attribute3,attribute4

what I want:

attribute1,attribute4
attribute2,attribute3,attribute4

can somebody help me ?


texts: Google Earth Engine: Name parameters in custom function

texts:
Situation
I'm working on a function with many parameters, and would like to be able to name them explicitly when calling the function.
I would usually call a function like this:
function get_sum(value_1, value2) {
var Sum = ee.Number(value_1).add(ee.Number(value2));
    print(Sum);
}
// Working
get_sum(75, 82);

Naming the parameters when calling the function leads to an error:
//Produces error
get_sum({value_1: 75, value2: 82});

Question
How do I do this? Is it something to do with served vs. client side?
Link to code


texts: QGIS Server GetCapabilities does not show map layers

texts: I install QGIS2.14 and QGIS Server on CentOS7.8. I make a map project with QGIS, then I put the qgs project into /var/www/wms, but when I access the GetCapabilities url through browser, but there is no layer info. What should I do?


texts:
I'm currently working on a forestry project for my thesis. It consists of the optimization of skids tracks planification using ALOS images.  To this end
I would like to extract obstacles such as steep slopes, stream and ridges as they are key zones for the skidding.

This question will focus on the extraction of ridges line.

I tried the inversion of the DEM and the hydrology analysis (Catchment Area then Channel Network using less-than 200 threshold methods) but the results aren't satisfying enough.

Do you know any way to extract ridges using QGIS and DEM?


texts: OpenLayers getText() from style function

texts:
I have a feature with a label set via feature.style.text. In my unit tests, I need to grab the value of the label.

Normally I would use getStyle().getText().getText(). However, in the example below I am using a style function, so calling getStyle() on my iconFeature returns a function and gives the error:

  iconFeature.getStyle(...).getText is not a function

var iconFeature = new Feature({
    geometry: new Point([250, 1000])
});

iconFeature.setStyle((f, r) => {
    return new Style({
        text: new Text({
            offsetX: 100 * (0.04 / r),
            padding: [0, 0, 0, 1000],
            text: 'Icon Label Vaue'
        })
    });
});

console.log(iconFeature.getStyle().getText().getText());

Any thoughts on how to get the label text of my feature's style?

In case it matters - This is an Angular app using OL 5


texts: ArcMap: Joining Tables for the Labeling Manager

texts:
I have two tables I would like to join, in an mdb, and have the label manager display columns from the joined table. I am able to do this through the ArcMap Gui (right clicking a layer and making a join through "Joins and Relates"), but I would like to do this programmatically. Is this possible?

I joined the two tables by following this example: http://resources.arcgis.com/en/help/arcobjects-net/conceptualhelp/index.html#//0001000002zw000000 but the label manager does not seem to recognize this join and as a result I am missing columns.

Thanks

Edit: I am working in C#.


texts: Creating OGC Services via GeoServer automatically?

texts:
I am new to GeoServer and wonder how I can create OGC Services via Geoserver automatically. I have to publish my data (Shapefiles) with a Java program as OGC Services. How can I do that?

I checked out the source code of GeoServer, but I don't know how to package as right jar file to use.


texts: Is it possible to create contour lines based on part of the points and not the entire dataset?


texts: Open Layers squished tiles until page zoomed

texts:
Upon initial load the map tiles are squished or sometimes do not appear properly as per image below. This is fixed when the page is zoomed in and out (not the map).

This is when the page is viewed on a 2K screen, any suggestions for how to fix this?


texts: TileMill-style based on the country the objects belong to

texts:
Is it possible to have different TileMill stylings in different countries? I need for instance different colours for highways based on the country the highway is in.
Is there a best practice solution to this problem? I know it's possible to store the borders of all countries in PostGIS an check the intersection of an object with all countries. I could also precalulate the corresponding country for all objects like in Nominatim. Are there any better solutions?


texts: Converting JSON structure gotten from Tracker to Esri JSON structure


texts: Joining attributes by location error when finding polygons within district boundaries?

texts:
I am using QGIS 3.5. I have a large shapefile (760 MB) with 1800 polygons representing different bird species ranges (download link)

I have a second shapefile with polygons of district boundaries (download link)

I have ensured the geometries are valid in both shapefiles and their CRS are equal. My goal is to identify the bird species present in each district. The output attribute table should include a column for the district ID (called c_code_11 in the data) of each district where the species exists.

I did a Join Attributes by Location with input layer as the species ranges and join layer as the district boundaries. I put 'intersects' and 'contains' as the geometric predicate, 'fields to add' as c_code_11 (the district ID), 'join type' as one-to-many, and then saved the joined layer as a shapefile to my local machine.

The output shapefile gets created on my computer, but at 3% completion, I get a list of errors saying the features could not be written to the shapefile. The process still runs to 100% completion (takes forever) and outputs a 4GB shapefile, but it is empty.

Any idea what is happening and how I can fix this?

Since the output shapefile is beyond the size limit, I re-ran the join attribute by location and saved the output as a GeoPackage. Now I am getting the following error.


texts: Saving vector layer with specific symbology using QGIS?

texts:
I am creating a polygon feature layer with soil classes in QGIS.

I need to share this layer with some colleagues and I need the layer to maintain the color ramp I have created.

Is there any specific format where the color ramp is saved together with the layer?

Since I don´t know how to do that, I was thinking on saving the color ramp that I created as a style.

How can I export this style?


texts: WMS with Access Control (ACL) for different users? Or solution for thousands of points on web map?

texts:
I have a set of sites in the thousands.  Not all users see same the same Sites on the map - thus ACL.  Data is stored in a non-GIS database (mssql, ora, etc) and right now is retrieved with a standard SQL query and returned to client (OpenLayers) with a JSP.

The current code is already optimized to do client-side clustering (so less of a mess on the client, but all the data is there), bounding-box retrieval (so a sub-set of data but have to go to server for every pan/zoom for new data) and even server-side clustering (to reduce data shipped to client).

What I would like to instead is generate images/tiles for all my data, similar to how Google Maps shows tons of little red dots for all the results - http://bit.ly/d73qrw [google maps search for "coffee"] and display on the client.  And when clicked-on do a quick ajax call to a WMS service to get info.  That's the idea at least.

But here's the problem - I can't just setup a standard WMS service in front of my data because not all users see same data.  Is there a way to generate these tiles on the fly or make WMS ACL-aware?

edit - 9/22/2010 - So I found out how google generates the tiles, or rather the technology behind it.  They are using Google Fusion Tables.  Store N rows in their tables and then if the data is lat/lon aware, the Fusion Table product can generate the tiles on the fly (!).  This is the kind of thing I'm looking for - performant map with tons of data.  But of course I need it to be ACL controlled.  Does it make sense to write a custom light-weight implementation of the WMS spec or modify an existing product?  Though geoserver seems like a lot to "just" modify to support ACL.

edit - 9/27/2010 - Some more info since adding bounty.  My data is in Oracle.  Ora spatial is not enabled.  Right now the data is extracted at the business level and converted to data, sent to client where the client puts the "dots" on the map.  ACL is done at the business logic level, not DB or ActivDir or anything like that.  Authentication is simple, but Authorization is not and so had to be captured in code.  Would like to know how to best create a WMS service to serve up 1000's of "dots" on the map where each user will see a different subset of dots.  Is the answer a CQL_FILTER?  But then how are the parameters set?  One idea that I have now is to do a 2-step process. First run in-house query to get list of ID's that user allowed to see, then build a WMS request string with those ID's in CQL_FILTER parameter.  Is there anything that simplifies this process?  And if I go with this, how can I add this layer as a "WMS" layer to an Open Layers client since to OL the end-point is my code to get IDs from DB not the actual WMS service on e.g. GeoServer?


texts: Arctoolbox Spatial join returns null value

texts: I am trying to join a polygon shapefine with DWG annotation, with ArcGIS Desktop 10.1 but all the fields returned are Null !!


texts:
I am running ArcGIS 10.2.2 Version 10.2.2.3552, using ArcMap.

I've watched the tutorial here (http://training.esri.com/Courses/ts_PythonAdd-Ins/player.cfm?) and I'm still not understanding this.  I am still in the beginner phase of understanding and using python.  I have created a script that when ran outside of arcmap everything works great.  Its complicated though to have people install pyscripter and then load my code.  My understanding is that i can create a button in arcmap that will run the code for them?  Is there a limit on what can be accomplished with this button?

When i post my script underneath and install via the add-in manager:

import arcpy
import pythonaddins

class ButtonClass4(object):
    """Implementation for Export Button_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):

I get a nice little missing icon with a circle and line going through it.  What am i missing here?


texts: Change total length unit using calculate geometry

texts: I want to change the total length from decimal degree to meter. How? the length was disable.


texts: Seeking tools for regionalization/redistricting?

texts:
I am looking for tools for regionalization / districting in GIS.

So far I have been able to identify

Districting ArcGIS extension (now also online)
REDCAP (currently limited to 50k input polygons)
iRedistrict (not yet released)
R method

Is there anything else available out there? Especially something that will have no problems in dealing with quite large dataset [~ 1.3 million polygons]?


texts: Using Open Universal Data tool correctly in MapInfo Pro?

texts:
After I lost settings file (MapInfow.prf), I can't properly open data using tool Open Universal Data. Which earlier worked without problems. Now when I try to open the dwg file, the layers are displayed incorrectly.

The layers are covered by the other, buildings are covered by parcels. I'm doing exactly what I always do. What should I setup in the options to correctly open dwg files.

I'm using MapInfo Pro 15.0 32-bit


texts:
Okay, so the Geoserver documentation (http://docs.geoserver.org/latest/en/user/configuration/status.html#config-serverstatus) says that clicking on the “Free Memory” button, cleans up memory marked for deletion by running the garbage collector.

So my question is why have this button in the first place?  Why isn't the GC being run automatically?  I am having to periodically click this button to improve performance, and that seems kind of lame.

I'm using the Concurrent Mark Sweep garbage collector as per the Geoserver "Container Considerations" help topic recommendation: http://docs.geoserver.org/latest/en/user/production/container.html#optimize-your-jvm