texts: QGIS 3 installation fails due to wrong python found

texts:
I try to install QGIS as seen in https://github.com/qgis/homebrew-qgisdev/issues/64 by executing:

brew install --no-sandbox qgis/qgisdev/qgis3-dev

but fail with:

CMake Error at /usr/local/Cellar/cmake/3.11.0/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find PythonInterp: Found unsuitable version "1.4", but required
  is at least "3" (found /usr/local/bin/python3)

How can I fix this, as

python --version
Python 2.7.14

geoheil@eduroam-078-104-029-050 ~                                                                                [17:34:02]
> $ python3 --version
Python 3.6.5

already returns the right versions of python?

full error

brew install --no-sandbox qgis/qgisdev/qgis3-dev
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 1 tap (homebrew/core).
No changes to formulae.

==> Installing qgis3-dev from qgis/qgisdev
==> Cloning https://github.com/qgis/QGIS.git
Updating /Users/geoheil/Library/Caches/Homebrew/qgis3-dev--git
==> Checking out branch master
==> /usr/local/bin/pip3 install --user future psycopg2 python-dateutil httplib2 pytz six nose2 pygments jinja2 pyyaml requests owslib
==> cmake -G Ninja -DCMAKE_C_FLAGS_RELEASE=-DNDEBUG -DCMAKE_CXX_FLAGS_RELEASE=-DNDEBUG -DCMAKE_INSTALL_PREFIX=/usr/local/Cellar/qgis3-d
Last 15 lines from /Users/geoheil/Library/Logs/Homebrew/qgis3-dev/02.cmake:
--   skipping QCA OpenSSL plugin C++ check
-- Found Libtasn1: /usr/local/include
-- Pedantic compiler settings enabled
-- Debug output enabled
CMake Error at /usr/local/Cellar/cmake/3.11.0/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find PythonInterp: Found unsuitable version "1.4", but required
  is at least "3" (found /usr/local/bin/python3)
Call Stack (most recent call first):
  /usr/local/Cellar/cmake/3.11.0/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:376 (_FPHSA_FAILURE_MESSAGE)
  /usr/local/Cellar/cmake/3.11.0/share/cmake/Modules/FindPythonInterp.cmake:155 (FIND_PACKAGE_HANDLE_STANDARD_ARGS)
  CMakeLists.txt:696 (FIND_PACKAGE)

-- Configuring incomplete, errors occurred!
See also "/tmp/qgis3-dev-20180411-72396-1w7hugc/build/CMakeFiles/CMakeOutput.log".

If reporting this issue please do so at (not Homebrew/brew or Homebrew/core):
https://github.com/qgis/homebrew-qgisdev/issues

These open issues may also help:
Update qgis3-dev.rb https://github.com/qgis/homebrew-qgisdev/pull/63
Update qgis3-dev.rb https://github.com/qgis/homebrew-qgisdev/pull/61
Python 3.6.0 upgrade breaks existing qgis3-dev https://github.com/qgis/homebrew-qgisdev/issues/11

Error: A newer Command Line Tools release is available.
Update them from Software Update in the App Store.


texts: Passing a stand alone script in a plugin

texts:
I have a few stand alone scripts and I am in the process of making a toolbar from them.

Most of my script have arguments like this :

##Points=vector
##Distance= number 2000
...

Those lines make a parameter window appear before the processing is done. This is the standard interface for stand-alone scripts.

Yet, is it possible to find a way to just link them in the plugins or do I have to make a new gui. Every button should be linked to one script.

Right now, the button have been done this way :

def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    self.add_action(
        ':/plugins/Concentration/param.png',
        text=self.tr(u'Param'),
        callback=self.run,
        parent=self.iface.mainWindow())

    self.add_action(
        ':/plugins/Concentration/group.png',
        text=self.tr(u'Algorithmes de groupement'),
        callback=self.run,
        parent=self.iface.mainWindow())
...

So if I could just find away aroud the callback=self.run to run a script just like I do in qgis(like callback = Param.py), it will save me the time to remake an interface for those scripts.


texts: How to calculate the area of a raster within a polygon which is greater than a specified value?

texts: I'm using QGIS (2.8.1-Wien) and I have a raster layer showing chlorophyll concentrations in the ocean. I have also defined number of regions of the ocean using polygons. What I'd like to calculate for each region is the area within that polygon where the raster is greater than X. For example, what area in each polygon has a chlorophyll concentration greater than 2. Any suggestions?


texts:
I have setup a consol app in C# that uses the GDAL.CORE NuGet,
it seems to work fine when converting most of the .shp files I have,
however with one .shp file I am getting the error below:

System.ApplicationException: 'EPSG PCS/GCS code 4326 not found in EPSG support files.
Is this a valid EPSG coordinate system?'

I am using the gcs.csv file below which seems to include this coordinate code = 4326

Also I have setup the path to my gdal-data correctly, so I know its reading my gcs.csv file.

Gdal.SetConfigOption("GDAL_DATA", "C:\\Program Files\\GDAL\\gdal-data");

Do you have any suggestions?


texts: Calculating distance of all points to central point using QGIS

texts:
I'm using QGIS and the points on this map represent different periods of the day in which coyotes were observed.  I am trying to figure out if coyotes are active at different times of the day based on their geographical distance from the center of the downtown area.
Is there a way to calculate the average distance of all the points in a shapefile from a central point?


texts: How can I create a form so the user has to be authenticated to access a project?

texts:
I have a login form in QGIS and when it is loaded, python gives me an error for, apparently, not recognizing the findchild part. This is due to the fact that the form isn't loaded before the script is executed, so there are no "childs". How can I connect the form so that the user has to be authenticated first in order to access the project?

import psycopg2
import hashlib
from PyQt4.QtCore import *
from PyQt4.QtGui import *

nameField = None
myDialog = None
passField = None

def formOpen(dialog,layerid,featureid):
global myDialog
myDialog = dialog

global nameField
nameField = dialog.findChild(QLineEdit,"username")

global passField
passField = dialog.findChild(QLineEdit,"password")

loginButton = dialog.findChild(QPushButton,"btnLogin")
loginButton.clicked.connect(authenticate)

def authenticate():

conn_string = "host='CORFOU' dbname='addisababa' user='addisababa' password=1"
conn = psycopg2.connect(conn_string)

username = nameField.getText()
hashpassword = hashlib.md5(passField.getText())

cursor = conn.cursor()
cursor.execute('SELECT * FROM users WHERE username=%s AND password=%s', (username, hashpassword))

if cursor.rowcount == 1:
myDialog.accept()
else:
msgBox = QMessageBox()
msgBox.setText("Wrong username or password")
msgBox.exec_()


texts: Is there any way of making label size dependant from the value of an attribute?

texts:
I know QGIS has a feature to make labels size dependant from an attribute:

The problem is that I don't know how to choose a maximum size or even a size variation-scale. It makes the size proportional to the values of the chosen attribute, which sometimes makes the lables too big (or too little):

I've tried to use the "expression" labeling (for example using logarithms), but I don't know how to make the maximum and minimum be less "extreme".

Is there a way (maybe through an expression) to make the size of the labels more visual-friendly? Something to make those labels "less proportional" to the attribute.


texts:
I want to create a group layer (or similar) of many raster files that only display when I tick them. They must not all be displaying when I make the layer, or my standard office entry level machine will crash.
I will need to recreate this "layer" as files are moved around and new ones are found on the network.
So I am after tricks for

a quick way to get them all into one "layer"
but with each raster being unticked/not showing so my machine does not crash

I have collected the list of 100 odd raster file paths in a text file. Can I easily turn this into some kind of text file that imports into ArcMap and it displays the list of grids without any being ticked to display? A text file with something wrapped around each file address, such as XML?
It does not matter if it converts into the ESRI binary format of a layer file. Making a large layer without dragging and dropping images is more important.
Added 10/9/2020
Follow Michael Stimson suggestion, I tried the following
import arcpy

mxd_path = arcpy.GetParameterAsText(0)
mxd = arcpy.mapping.MapDocument(mxd_path)

text = arcpy.GetParameterAsText(1)

layer_path = arcpy.GetParameterAsText(2)
lyr = arcpy.mapping.Layer(layer_path)

df = arcpy.mapping.ListDataFrames(mxd)

log = open(text + ".log", "w")

f = open(text, "r").readlines()
for filepath in f:
    log.write(filepath)
    newlayer = arcpy.mapping.Layer(filepath)
    arcpy.mapping.AddLayerToGroup(df, lyr, newlayer)

lyr.save()

The only failure is AddLayerToGroup. The UNC network address successfully creates a layer and I even isolated to check the lyr.save line works by itself.
Error:
Traceback (most recent call last):
  File "K:\Tim K\Program Resources\ArcGIS\Arcgis10.1\GridsToLayer.py", line 22, in <module>
    arcpy.mapping.AddLayerToGroup(df, lyr, newlayer)
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\mapping.py", line 87, in AddLayerToGroup
    assert isinstance(data_frame, DataFrame)
AssertionError

Failed to execute (GridsToLayer).

I suspect the problem is in the AddLayerToGroup help where it states the layer I want to add to the group layer must be "a layer file on disk or a layer file in a map document". That would mean copying the raster files to new disk location and converting each into layer. I was hoping to simply reference where they are now. That location may change a bit each week a few files pass through a QA/QC "process" by others.
It looks like I want an unmanaged Raster Catalog, which references files instead of copying to a gdb. But, can it reference a workspace for every raster?
Turns out yes, but its not displaying the rasters, just polygons where they are located.
See answer below - changed to making an mxd of rasters.


texts:
I have a function I am creating for some iterative task I need to do but I cannot get the geometry coming out of a record data type into the execute properly:

DO
$$
DECLARE
records record;
images record;
BEGIN

FOR records IN SELECT opening, shape FROM d.blocks_2013_2018
LOOP
    TRUNCATE d.xx_delete;

    INSERT INTO d.xx_delete (opening_number, image_name)
    SELECT records.opening, a.image_name
    FROM d.image_tiles_with_blocks a
    WHERE ST_Intersects(a.shape, records.shape);

    FOR images IN SELECT * FROM d.xx_delete
    LOOP
        EXECUTE 'INSERT INTO d.cc_analysis (opening_number, t_id, ht, d_huang, d_huang_adj, d_cortini, geom)
                 SELECT '''||images.opening_number||''', b.t_id, b.max, b.d_huang, b.d_huang_adj, b.d_cortini, b.geom
                 FROM p839.'||images.image_name||'_t b
                 WHERE ST_Intersects(%, b.shape)', records.shape;
        COMMIT;
    END LOOP;
END LOOP;
END;
$$
;

It is the record.shape at the end of the code that I cannot get into the code properly. Currently, I am getting an error stating that the "query ... returned two columns"


texts: Examples of web based Overweight/Oversize load permit routing on County Roads

texts:
Does anyone know of any US counties that have set up a web based overweight/oversize load permit routing system, preferably one that leverages GIS?

I know many state DOT's are doing this, but I don't see any that integrate with county maintained roads.


texts: How to get all fields of selected features in C#

texts:
I am looking for a way to get all selected features from a map which belongs to a specific feature class (shapefile, geo database,…).

To be more specific and to give you a real case for my question: I want check if the feature class of selected feature has a special field. If this is the case I want to get the value of this field from the selected feature and add this value into a combobox. This process should be repeated for all selected features.

I tried the following but it crashed with an object reference error.

List<string> list = new List<string>();
IMap Pmap = doc.FocusMap;
IEnumFeature pEnumFeat = (IEnumFeature)Pmap.FeatureSelection;
pEnumFeat.Reset();
IFields fields;
try
{
  IFeature pfeat = pEnumFeat.Next();
  while (pfeat != null)
  {
    fields = pfeat.Fields;
    int x = fields.FindField("ID_K");
    if (!x.Equals(-1))
    {
      list.Add(pfeat.get_Value(x).ToString());
    }
    pfeat = pEnumFeat.Next();
    }
}
catch (Exception e)
{
  Console.WriteLine(e.Message);
}

UPDATE:

I've found a solution. Here is the working code:

List<string> list = new List<string>();
IMap Pmap = doc.FocusMap;
IEnumFeature pEnumFeat = (IEnumFeature)Pmap.FeatureSelection;
IFields fields;
pEnumFeat.Reset();
IEnumFeatureSetup enumFeatSetup = (IEnumFeatureSetup)pEnumFeat;
enumFeatSetup.AllFields = true;
try
{
   IFeature pfeat = pEnumFeat.Next();

   while (pfeat != null)
   {
      fields = pfeat.Fields;
      int x = fields.FindField("ID_K");
      list.Add(pfeat.get_Value(x).ToString());

      pfeat = pEnumFeat.Next();
    }
 }
 catch (Exception e)
 {
    Console.WriteLine(e.Message);
 }

These two line are especially important.

IEnumFeatureSetup enumFeatSetup = (IEnumFeatureSetup)pEnumFeat;
enumFeatSetup.AllFields = true;

If you don't use them you've only the shapefile included (default) and not the other fields and values.


texts: ModuleNotFoundError when installing plugin from zip

texts:
I have a plugin that I have developed for work and we have less than 10 users that will want to use it, so rather than publishing on the QGIS online plugin repository (https://plugins.qgis.org/plugins/), I thought I could zip the plugin directory and disseminate that instead, making use of the 'INSTALL from ZIP' option on the plugin manager in QGIS.

I thought it would be as simple as zipping the directory containing the plugin code and simply pointing to that file, but I'm getting a 'ModuleNotFoundError':

Couldn't load plugin 'datacube_qgis_plugin/.git'

ModuleNotFoundError: No module named 'datacube_qgis_plugin/'
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.8/apps/qgis/./python\qgis\utils.py", line 312, in loadPlugin
    __import__(packageName)
  File "C:/PROGRA~1/QGIS3~1.8/apps/qgis/./python\qgis\utils.py", line 737, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'datacube_qgis_plugin/'

My directory structure looks like:

datacube_qgis_plugin.zip
                   |--- > datacube_qgis_plugin/
                                              |---> .py files etc.

Is there an exact method for 'zipping' a plugin for QGIS import?

ADDITIONAL INFO - If I exit the box displaying the error above, it then proceeds to say 'plugin installed successfully' - so why the error?


texts: How do I get a .csv file of thousands of GPS coordinates along a stretch of a highway? (GPS Trace?)

texts: I need a large amount of [latitude, longitude] data points along Texas' I-35 from Hillsboro to Salado. Essentially need a GPS Trace with minimum ~1,000 points.


texts: FMRC feature in THREDDS presents timesteps incorrectly

texts:
I'm trying to set up a Forecast Model Run Collection using THREDDS. I have a collection of netCDF files, where the time at which the simulation was run is part of the filename:

forecast_YYMMDDHHmm.nc

In the files, the time variable is an integer, and gives time as number of hours since a given date. This date is the same in all my files. Since the times are given in integer hours, there should be absolutely no uncertainty as to whether two timesteps are the same.

However, when the THREDDS server presents this as an FMRC using the "best" model (for a given timpstep, it always serves data from the most recently performed simulation), it converts the time variable to a double, and somehow manages to present some timesteps which are actually the same as slightly different. For example, all of these are actually the same in my data, yet are presented as three different timesteps:

-65.10000033333334, -65.10000033333333, -65.1

The reason that the numbers are not whole hours is that the server also changes the reference point for the time axis from that specified in my files, to the simulation run time of the first file. This is fair enough, if only the data would be presented correctly.

I find the THREDDS documentation to be a collection of nice examples rather than an exhaustive list of all requirements and all available options, so I haven't been able to figure out how to fix this problem.


texts: QGIS Offline Editing Plugin: Automate manual steps to sync offline database with cloud database

texts:
The Offline Editing Plugin offers a good way to work offline on a project and sync it later with a cloud database. The problem is that before I go offline I need to explicitly do some manual steps to convert the project to an offline project. That means if I'm already offline I can't start offline editing.
How could I solve this problem? My first idea was to automate this process and whenever I sync the project back with the cloud database I directly convert it back to an offline project.
Even better would be if QGIS detects if I have a connection to the cloud. And if so I can work directly on the cloud database. And if I have no connection It will save it in the offline project an as soon as it has a connection to the cloud it will sync it automatically.
Does anyone have experiences with a similar use case? Or any advice on how to solve this problem the best?


texts: Can I export a list of user and their roles from ArcGIS Server 10.0?

texts: I need to export a list of all users and roles in ArcGIS Server 10.0. Does anyone know if there is a method to do this?


texts:
I try to use this code from QGIS 2.x version in QGIS 3.0:

li = iface.legendInterface()
li.setLayerVisible(myLayer, False)

I found QgsAnnotation but can not use it.


texts:
I just extracted some river paths from some old maps I acquired a few weeks ago, they seem accurate enough and will be very interesting to put against my theoretical mapping calculated by using the mnt.

The little issue I have is as follows: Even if I am satisfied by what was extracted and turned into a linestring by SAGA, a lot of gaps are still present. This is due to the fact that some bridges and roads overlap with the rivers and the accuracy of my very own pixel extraction. I am looking for the best way to fill those gaps, using Postgres if possible.

I tried a few things already but it wasn't very useful or accurate. I tried using centroids, startpoints, endpoints combined with Delaunay triangulation and Voronoi polygons. It did simplify the geometry and filled some of the gaps but the query was very complex and convoluted, even after optimization it took a lot of time to do anything even on only a few segments. That and I don't think I am going about it the right way.


texts: Use a PNG image as a correctly georeferenced background in order to trace lines from it?

texts:
I have a PNG image which looks like this and has some lines which I need to trace (the green and red lines):

The coordinate system of the map from which the image is taken is EPSG:3034. How can I position the image correctly on a shapefile such as this one (a shapefile with coordinate system EPSG:32632) using QGIS:

I tried to georeference the PNG file (with the Georeferencer plugin) by defining the coordinates of four points, but somehow I am getting a distorted final image. As an example, when I am inputting the East and North coordinates in the georeferencer (QGIS) I use the E/N coordinates of the form 44.214555, 22.67459 (WGS84).

How do I solve the problem?

The idea is to be able to extract the information regarding the location and length of the lines from the PNG file in a new EPSG:32632 shapefile.

The result I get, when georeferencing 4 points (geographical extremities of the country. i.e. most northern point, most southern point, etc.):


texts: Creating new layer in ArcGIS for Desktop with Python Toolbox?

texts:
I would like to create new shapefile layer or memory layer with one point with certain coordinates in ArcGIS 10.2 with Python Toolbox. Which function should I use? I wrote this code:

import arcpy

class Toolbox(object):
    def __init__(self):
        self.label = "Python Toolbox"
        self.alias = ""
        self.tools = [XY]

class XY(object):
    def __init__(self):
        self.label = "XY"
        self.description = "Coordinates to point"
        self.canRunInBackground = False

    def getParameterInfo(self):
        parameters = []

        x = arcpy.Parameter(
            displayName="X",
            name="x",
            datatype="String",
            parameterType="Required",
            direction="Input")
        parameters.append(x)

        y = arcpy.Parameter(
            displayName="Y",
            name="y",
            datatype="String",
            parameterType="Required",
            direction="Input")
        parameters.append(y)

        return parameters

    def isLicensed(self):
        return True

    def updateParameters(self, parameters):
        return

    def updateMessages(self, parameters):
        return

    def execute(self, parameters, messages):
        x = float(parameters[0].valueAsText)
        y = float(parameters[1].valueAsText)


texts:
I am attempting to map the flow of remittances into Japan using the FlowMapper plugin in QGIS. When I use the Style--> Advanced--> Size Scale field--> magnitude command, the lines become so thick and varied that they are impossible to distinguish from one another (in fact my map becomes one giant block of "lines"). Does anyone have any advice on how to fix this other than manually changing the size of lines classified under equal intervals?

Is there a way to label the flow lines by country of origin so that the text is in and oriented with the flow line, not located on top (that is, once I'm able to get the flow lines to their desired thickness)?

How do I make the lines curved so they are more aesthetically pleasing?


texts: gdal.py as part of bat script in windows 7

texts:
I need to run the following bat file from a MS4W shell which has gdal_retile.py as the final option. I have added C:\ms4w\python\gdal to my pythonpath but within a pyscripter shell I can not import gdal or gdal_retile.py

Without a linux box is there a way to get the final step to work?

CODE of P.bat File

mkdir -p %2\%3
gdal_translate -expand rgba -of GTiff -a_srs EPSG:3857 -co TILED=YES -co BIGTIFF=YES %1 %2\%3\%3_translated.tif
gdalwarp -multi --config GDAL_CACHEMAX 500 -wm 500 -r near -co "COMPRESS=DEFLATE" -co "TILED=YES" -co "BIGTIFF=YES" -srcnodata 255 -dstnodata 255 -overwrite %2\%3\%3_translated.tif %2\%3\%3_rebuilt.tif
mkdir -p %2\%3\pyramids
gdal_retile.py -v -r bilinear -levels 4 -ps 2048 2048 -co "TILED=YES" -co "COMPRESS=DEFLATE" -targetDir %2\pyramids %2\%3\%3_rebuilt.tif

Pause

OUTPUT of process

  GDAL, mapserv, mapcache, PROJ, and shapelib dll paths set GDAL_DATA
  path set GDAL_DRIVER_PATH set PROJ_LIB set CURL_CA_BUNDLE set
  C:\Users\georgec\Desktop>P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer
  _Creation\4Band_translated_pyramid\p.bat P:\2012\183_TownPlanning_Symbology\Outp
  ut\8_Bit_TIF_Colourmap_RS\TRC_WC_WaterResource_Catchment\TRC_WC_WaterResourceCat
  chment.tif
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\4Ban
  d_translated_pyramid wc

  C:\Users\georgec\Desktop>mkdir -p
  P:\2012\183_TownPlanning_Symbology\Working\Ras
  ter_Layer_Creation\4Band_translated_pyramid\wc

  C:\Users\georgec\Desktop>gdal_translate -expand rgba -of GTiff -a_srs
  EPSG:3857
  -co TILED=YES -co BIGTIFF=YES P:\2012\183_TownPlanning_Symbology\Output\8_Bit_TI
  F_Colourmap_RS\TRC_WC_WaterResource_Catchment\TRC_WC_WaterResourceCatchment.tif
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\4Band_translate
  d_pyramid\wc\wc_translated.tif Input file size is 19844, 26300
  0...10...20...30...40...50...60...70...80...90...100 - done.

  C:\Users\georgec\Desktop>gdalwarp -multi --config GDAL_CACHEMAX 500
  -wm 500 -r n ear -co "COMPRESS=DEFLATE" -co "TILED=YES" -co "BIGTIFF=YES" -srcnodata 255 -dst nodata 255 -overwrite
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Cr
  eation\4Band_translated_pyramid\wc\wc_translated.tif
  P:\2012\183_TownPlanning_Sy
  mbology\Working\Raster_Layer_Creation\4Band_translated_pyramid\wc\wc_rebuilt.tif

  Creating output file that is 19844P x 26300L. Processing input file
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Cr
  eation\4Band_translated_pyramid\wc\wc_translated.tif. Using band 4 of
  source image as alpha. Using band 4 of destination image as alpha.
  0...10...20...30...40.ERROR 2: Out of memory allocating 130520448 bytes for Unif iedSrcDensity mask. ..50..
  C:\Users\georgec\Desktop>mkdir -p
  P:\2012\183_TownPlanning_Symbology\Working\Ras
  ter_Layer_Creation\4Band_translated_pyramid\wc\pyramids A subdirectory
  or file -p already exists. Error occurred while processing: -p.

  C:\Users\georgec\Desktop>gdal_retile.py -v -r bilinear -levels 4 -ps
  2048 2048 - co "TILED=YES" -co "COMPRESS=DEFLATE" -targetDir
  P:\2012\183_TownPlanning_Symbol
  ogy\Working\Raster_Layer_Creation\4Band_translated_pyramid\pyramids
  P:\2012\183_
  TownPlanning_Symbology\Working\Raster_Layer_Creation\4Band_translated_pyramid\wc
  \wc_rebuilt.tif Traceback (most recent call last):   File
  "C:\ms4w\python\gdal\gdal_retile.py", line 38, in
      import gdal ImportError: No module named gdal

  C:\Users\georgec\Desktop>Pause Press any key to continue . . .

I also get the following Error

  ERROR 2: Out of memory allocating 130520448 bytes for Unif iedSrcDensity mask


texts: Linking census shapefiles with geocoded addresses using QGIS

texts: I'm trying to merge a census shapefile of New Orleans with a health outcomes dataset by address with QGIS. I have geocoded the addresses and added them as a layer in QGIS but don't know on which variable within the street level city shapefile to merge them.  The dataset has long/lat but the shapefile does not.


texts: ArcPy rollback update / insert cursors if error occurs

texts:
Is it possible to rollback an update made with an update cursor or an insert made with an insert cursor if an error occurs?

I need to go through a feature class and break numbers up into 100 groups. e.g.

num_start | num_end | other_value
---------------------------------
100       | 300     | foo

Needs to be:

num_start | num_end | other_value
---------------------------------
100       | 199     | foo
200       | 299     | foo
300       | 300     | foo

My update cursor changes the first row, and then the insert cursors inserts the remaining rows. The update and insert work fine, but when an error occurs the updated values remain causing the sequences to be incorrect, which of course means I need to re-export the feature class and start over again.

I am using Python 2.7.13.

Libraries:  These are the libraries I import

import arcpy
from arcpy import env
import sys
import os
import logging
import time
from functools import partial # Allows function calls with parameters inside of dictionaries

What I tried: I though I could have a variable set before a try statement and then reset the data in an exception.

with arcpy.da.UpdateCursor(fc, "*") as updateCurs:
    oldValue = [] # initialize empty value
    for row in updateCurs:
        oldValue = [rowVal for rowVal in row] # read row into array
       try:
           # do stuff
           del updateCurs
       except Exception, e:
              # reset row to old row value
              for i in range(0, len(oldValue)):
                  row[i] = oldValue[i] # set row values
              # rollback
              updateCurs.updateRow(oldValue)

The problem here is that it would only correct the last row updated, but not the all the previous ones before it.

Row attributes: all other row attributes like geometry are set using the initial row being expanded, which currently are all points.


texts:
I have asked this question a while back but within arcgis environment without using addin. Now i am trying to do the same task to get pixel value when user click anywhere on screen, to use in later steps. Following pythonaddin button code snippet returns nothing. Any suggestions where is the catch?

class Val(object):
"""Implementation for Pythonaddin_addin.genW (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    geog = arcpy.GetParameterAsText(0)
    ras = arcpy.PointToRaster_conversion(geog, "ID", "in_memory/ras", "MAXIMUM", "",30)
    pass


texts: QGIS Modeler using File/Folder input paramerter

texts:
I would like to use the File/Folder input parameter to set a path for all output layer files. I was trying to connect the File/Folder parameter to my output from a model and set a default folder path.
I was trying to call the path by using @output_Directory but the path is not used. How can I connect the File/Folder input parameter with the model output?


texts: Using arcpy.da.InsertCursor to insert entire row that is fetched from search cursor?

texts:
This is a simple process using the legacy cursors, but I cannot figure out how to do it with the newer Insert Cursor from the Data Access module. I basically want to take the entire row from a Search Cursor created on an SDE feature class in one database, and insert that row into an SDE feature class in another database using the Insert Cursor. The code below is how I am doing this using the legacy cursors, and it works quite well, but I would like to take advantage of the faster performance of the Data Access cursors. This runs on a list of feature classes that all have different numbers of fields (and different geometry types), though the schema for each feature class is identical between the databases (the databases are basically copies of one another):

    sourceAddRows = arcpy.SearchCursor(sourceFC)
    targetAddRows = arcpy.InsertCursor(targetFC)
    for sourceAddRow in sourceAddRows:
        targetAddRows.insertRow(sourceAddRow)


texts: Extracting information from Google Earth Engine image by row and column position?

texts:
I would like to iterate over an image in Google Earth Engine, but not sure how to reference the location in an image by row and column. Here is some example code:

var ROI = ee.Geometry.Rectangle(-82.5574, 35.5958, -82.5563, 35.59515); // define area

var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterDate('2014-05-01', '2015-10-01')
    .filterBounds(ROI)
    .limit(1);

var composite = ee.Algorithms.Landsat.simpleComposite(collection);
var image = composite.clip(ROI);

Map.centerObject(ROI);
Map.addLayer(composite, {bands: 'B4,B3,B2', min:0, max:100}, "composite", false);
Map.addLayer(image, {bands: 'B4,B3,B2', min:0, max:100},"image");

// Extract projection information from the first band of a sample image in the collection.
var sample = ee.Image(collection.first());
var projection = sample.select('B1').projection();

// Reproject the composite image to the selected projection.
image = composite.reproject(projection).clip(ROI);

// Return the image metadata to the client.
var image_info = image.getInfo();

// On the clientside, using Javascript, print out the dimensions.
print('band B1 dimensions: ' + image_info.bands[0].dimensions);

print('image information', image);

The resulting image is 3 by 4. How can I reference the values at, say, position [0,0] or [1,3]?

I have tried using regular javascript syntax image[0][0], to no avail.


texts: Missing options on database connection in ArcCatalog

texts:
I have 2 databases on the same server. I can add connections to them both without a hitch in ArcCatalog (or catalog in ArcMap). The problem is that I seem to be losing some options when right clicking one of the database connections. Here is the first connection:

and here is the second:

As you can see, on my second connection, I do not have the option for "new", "import", or "export". Another user with seemingly identical permissions still has the options when they right click, so it only seems to be affecting me and only with this one database connection. How do I regain the ability to select these options? Also, I am on 10.1 for everything.


texts: Splitting polygon into two equal parts from specific point to polygon dividing line?

texts:
I want to split a polygon into two equal parts but I want from a specific point to the polygon dividing line.

How to do this using ArcGIS Desktop 10?


texts: Are there different definitions of "merge"?

texts:
The merge tool in ArcGIS ArcToolbox is used to combine two layers, but it seems that there is also a Merge option under Modify Features that is used to combine and dissolve adjacent features within a layer (https://pro.arcgis.com/en/pro-app/help/editing/merge-features-into-one-feature.htm).

Clearly the term is used at least two different ways, but are there two different technical definitions for "merge" within the context of GIS, or is one of these a less formal use of the term?


texts:
I have an image service connection which I'm trying to download individual tiles from.

Using ArcMap 10.5, I've selected a tile from the attribute table, right-clicked the image service > Data > Download Selected Rasters... and I get the following error message (see image below).

  Failed to retrieve the download files for the list request.

Am I doing something wrong related to the image service permissions?

The Esri help site does not address this error and I could not find a relevant forums that address this issue as well.


texts: Georeferencing non-projected pdf using bounding box in FME?

texts:
I have two similar looking and extend map files. Pdf -map (N4443) with vector but not projected and the raster tiff map (N4443) with projection.

I need to georeference pdf-vector objects using tiff -file bounding box.

I have non-georeferenced pdf -map (in vector format) and georeferenced raster map with bounding box. They have same bounding box but pdf- bounding box is not projected. How I can attach raster bounding box to using pdf also? What tools can be found in FME to help me?

error: 2019-03-22 14:51:11|   1.5|  0.0|WARN  |AffineWarper(WarpFactory): Must have at least 4 control points to perform an AFFINE transformation


texts:
I am trying to disaggregate an existing raster with a cell size of 10000 to a new one of 2000. The larger raster represents the UK in 10km squares, and I want the same extent but in 2km squares. I am using the disaggregate function in the raster package, but what should be simple task returns an error message Error in .local(x, values) : values must be a vector

The larger raster does according to the raster properties in ArcGIS have nodata values set to -32768, and I have tried changing these to 0

grid10[grid10==-32768] <- 0

# and in case there are NA's

grid10[is.na(grid10[])] <- 0

# but these make no difference

grid2 <- disaggregate(grid10, fact=5) # i.e 10000 to 2000 is 5 times smaller

  Error in .local(x, values) : values must be a vector

What I am I doing wrong?


texts: SQL Server 2008 Distance and Performance Problem

texts:
I had been working on PostGIS to do some geographic processes and due to some reasons, I had to start to use SQL Server 2008.
I ported same code for SQL Server but performance is almost dying when I sort according to distances (ORDER BY dist). I dont know if i am missing something on SQL server because i am new on it.

DECLARE @myPoint geography;
SET @myPoint = geography::STGeomFromText('POINT(28.8801844444445 41.0299561111109)',4326);
SELECT TOP 8 id,
@myPoint.STDistance(geom.ToString()) as dist
FROM MY_DB
ORDER BY dist

This query takes almost 4 seconds on SQL Server while it is taking 120 ms on postgis.


texts: Why is ArcPy Describe Extent overshooting the extents?

texts:
Why is arcpy.Describe(file.shp).extent overshooting the extent? Screenshot below. I have confirmed there are no tiny lines lingering out at the edges. I'm in UTM coordinates.

extent = arcpy.Describe(Lines.shp).extent
xmin = extent.XMin
xmax = extent.XMax
ymin = extent.YMin
ymax = extent.YMax

print xmin, xmax, ymin, ymax

Print output:

540375.245005 546724.102422 3814174.9483 3818495.05067

EDIT: I tried adding in a minimum bounding geometry, still same issue.

test_out = outpath + '\\' + 'temp_convexHull.shp'
arcpy.MinimumBoundingGeometry_management(Lines.shp, test_out, "CONVEX_HULL", "ALL", "", "")

extent = arcpy.Describe(test_out).extent
xmin = extent.XMin
xmax = extent.XMax
ymin = extent.YMin
ymax = extent.YMax

print xmin, xmax, ymin, ymax


texts: Using the Field Calculator in a Domain field

texts:
I have a coded domain field with 10 dropdown options.  Possible choices include integers and text.

I need to change a selection of 250 features (which have differing values) to the same value.  This value is one of the dropdown values.

When I use the field calculator (field="value") I get an error.  The calculated value is invalid.....

Anyone know how to change multiple values within a coded domain field?


texts:
I have a shapefile which contains a 20 hectare polygon. I want to demarcate 9,48 hectares from the entire polygon with a parallel.

I can make the parallel line with the longer line but I can not do the incision...

Can you help me how to do the incision?


texts: Clipping layer with every feature in another layer using ModelBuilder

texts:
Currently working with ArcMap 10.5, from a gdb, and using ModelBuilder. My python knowledge is null.

I have a road layer (line), and then a city limits layer (polygon) with various feauters (up to 30) each being a department of it. What I'm trying to achieve seems pretty simple but I am unable to find the way.

I want to clip the road layer with the city layer in a way that it creates a new road layer in the shape of every department. Something like this:

Road_1
Road_2
Road_3
Road_...
Road_29
Road_30

The next step would be to use Pivot Table on each one of these layers to know how many metres there are for each road level (1-5).

Is there any tool that does this kind of clip in a single process? I'm working on a ModelBuilder so I'd rather not have to add the same process 30 times unless it's absolutely necessary.


texts: Openlayers 4 - Clustering view - Select and Modify Interactions

texts:
I'm working on a specialized browser-based editing app based on openlayers 4. I'm displaying features (pictures, data Elements, docs) from a few source Layers which are geo referenced, using drag & drop, etc. and a server based datastorage.
Some of thes features are in the same location, others close together. These items should be selected and dragged individually, but while not selected, they should be display as a 'cluster' with e.g. a folder icon.

I'd like to have a hover which explode the cluster to display the original Features and a linked line to the original Point. And make the original features selectable.
Then, if an item is selected, the others should collapse to the cluster/folder-icon and the selected item shows on top of it.
If the item is deselected, then it disappears in the folder.
If the original Feature is moved or edited the original Data should be updated.

I took a look at the cluster source but ran in many difficulties. They first was: I'd like to cluster over a view vector Sources, with different behaviour and styles of their features. It wasn't that complicated to extend the ol.source.Cluster with the ability to handle a collection of VectorSources, including a addSource and a dropSource function (I'll make it available when it is finished and maybe useful), but I'm not happy with it. The double layer of features makes it confusing and creates the need to distinguish between clustered and unclustered features, nearly everywhere, resulting logicaly in an ol.clusterFeature object

I also created Select and Modify Interactions but they work only with plain VectorSources, but nut with clustered ones. Because there is the need to handle both, the clusterFeature and the originalFeatures. The cluster breaks the other mechanism of editing - or show me where I'm wrong...

So I'm not happy with this approach of handling displaying of closed located features.
While thinking of it, the thought grows to me that it is not a data source problem, but something related to the view and rendering.
Together with the question how to arrange the overlapping feature-display-areas.

Where is the Point to start in ol4, when I'd like to implement some clustering and arrange display-areas for viewing the data the described way?

Or what do you think is best practise to reach my goals with ol4?

Ervin Peters


texts: ERROR 000338 when appending two datasets into one with arcpy

texts:
In an enterprise geodatabase with three datasets I'm trying to append two of the datasets into the third, however I get an error which states that Inputs must be either all Feature Classes, Tables or Rasters; not mixed.. This seems odd as the datasets only consists of Feature classes.

The feature classes is structured with the same name but a different code suffix:
eg:
lsat_6865 : lsat_6866
lsar_6865 : lsar_6866
lsar_6865_anno : lsar_6866_anno

The number codes should be stripped from the name in the resulting dataset.

The entire error message below:

Traceback (most recent call last):
  File "C:\scripts\test\append.py", line 26, in <module>
    arcpy.Append_management([WorkspaceES + '/' + fc + MCodeES, WorkspaceGT + '/' + fc + MCodeGT], fc, 'TEST', '', '')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 4271, in Append
    raise e
ExecuteError:
ERROR 000338: Inputs must be either all Feature Classes, Tables or Rasters; not mixed.

import arcpy

MCodeES = '6865'
MCodeGT  = '6866'

WorkspaceES   = 'C:/scripts/sde_mcde.sde/sde_mcde.sde.ADO.ES_tab'
WorkspaceGT   = 'C:/scripts/sde_mcde.sde/sde_mcde.sde.ADO.GT_tab'

arcpy.env.workspace = 'C:/scripts/sde_mcde.sde/sde_mcde.sde.ADO.GT_tab.ADO.ES_GT.FK'

# Get a list of data sets and feature class
fcList = arcpy.ListFeatureClasses()
fcList.sort()

# Loop through the list
for fc in fcList:
    arcpy.TruncateTable_management(fc)
    if fc.partition('_anno')[1].strip() == '_anno':
        arcpy.Append_management([WorkspaceES + '/' + fc.partition('_anno')[0].strip() + '_' + MCodeES + '_anno', WorkspaceGT + '/' + fc.partition('_anno')[0].strip() + '_' + MCodeGT + '_anno'], fc, 'NO_TEST', '', '')
        print(fc)
    else:
        arcpy.Append_management([WorkspaceES + '/' + fc + MCodeES, WorkspaceGT + '/' + fc + MCodeGT], fc, 'TEST', '', '')
        print(fc)

print('/nPress Enter to continue ...')
raw_input()


texts:
I want to create a data frame in R combining the bands wth the shapefile. The final file will have the pixels as rows, the bands as well as the other columns from the attribute table from the shapefile.

I want a csv table where the columns are the Landsat bands (say band 2, 3 and 4) along with other columns that are there in the attribute table of the shapefile (Country, State, Village name, Lat, Long etc)


texts: problem in merging adjacent address polygons in qgis

texts: I am having 15 polygon shape files of a sub locality of a city that is having address in their attribute table. I want to merge these polygons based on the address in the attribute table?how can i do this in Qgis. I also attempted vector...>>>data management tools...>>>>merge shape files to one.But it is giving me the combined polygon but when I look into the attribute    table the values shows me null? What is the problem behind it?


texts: Activating custom ITool from form without adding it to ArcMap AddIn toolbar?

texts:
I'm working on an add-in for ArcMap 10.0 that adds a toolbar to ArcMap. One command (OpenModelessDialogCommand) button on that toolbar opens a modeless WinForms dialog, from which a tool (MyTool) can be activated in order to e.g. select a feature on the map.

My Config.esriaddinx contains these command and toolbar declarations:

<Commands>
  <!-- this is the command that opens the modeless WinForms form, from where
       MyTool is available: -->
  <Button id="OpenModelessFormCommand" ... />
  <!-- MyTool is not directly referenced in any toolbar defined in this file: -->
  <Tool id="MyTool" class="MyTool" ... />
</Commands>
<Toolbars>
  <Toolbar ...>
    <Items>
      <Button refID="OpenModelessFormCommand" />
    </Items>
  </Toolbar>
</Toolbars>

What I'm having problems with is activating MyTool in the form. All I've found on the internet is code samples along the lines of:

// get a reference to an instance of MyTool:
ICommandItem myTool = ArcMap.Application.Document.CommandBars.Find("MyTool");
 // activate MyTool:
ArcMap.Application.CurrentTool = myTool;

However, this apparently requires that MyTool actually appears in a command bar (e.g. toolbar) of my add-in. But that's not the case. So, I've tried this next:

ITool myTool = new MyTool();
ArcMap.Application.CurrentTool = myTool;  // Type mismatch! An ICommandItem is expected.

I've even looked into adding an invisible AxToolbarControl to my form and adding a button for MyTool there; but then I'm running into problems on how to connect that toolbar (via SetBuddyControl) to the opened document's map. I don't want the tool to function in a separate AxMapControl, I want it to work directly with the main map shown in ArcMap.

Question:
How do I activate a custom tool that is not added to any toolbar (or other command bar, for that matter)?


texts: Viewshed Analysis incorporating tree height

texts:
I'm trying to determine the visibility of a property (polygon) from surrounding sites within the same county.  Sounds simple, but here's the catch: it's a solar farm analysis project, and the county requires the subject site be 'invisible' or not visible from any other property in the county.  That said, the property currently has a thick loblolly pine plantation along its edges (trees are approximately 20m in height). So I've run a bare earth viewshed analysis (no problem with that one) BUT am trying to figure out how to run a separate analysis that incorporates/considers the height of the pines as an obscuring factor. Any chance you could help with this?

I don't have a DSM to work with, just a DEM.


texts: Why QGIS3 is not running processing models and scripts in background?

texts:
I have a problem with QGIS 3.6.0 on macos and QGIS3.2.3 on Windows7 with running processing scripts and models in background. According to QGIS3 docs they are run in background thread by default (if there is no FlagNoThreading flag defined.

Unfortunately there is no "run in background" button on processing dialog (just "Run") and they are run in main GUI thread (or in other thread but blocking GUI anyway).

I couldn't find any QGIS configuration setting for this. What could be the reason for QGIS3 to run them in GUI-blocking mode?


texts: How to run shp2pgsql in python without reporting every row loaded?

texts:
I'm running shp2pgsql in python with a simple script, which works just fine:

import os
import subprocess

cmd = 'shp2pgsql -s 4326 ../data/shp/statistical_neighborhoods.shp temp_table | psql -h hostname -d databasename -U username'

subprocess.call(cmd, shell=True)

The problem is, for every record loaded, there is a message reported:

INSERT 0 1

That message is reported for the 78 records in the neighborhoods layer in the example above, but for a file with 280,000 records (address point file from the City, or Parcels layer), that would be a very unnecessary amount of messages.

Is this problem coming from shp2pgsql, or is it a result of the subprocess.call() function?


texts: ogr2ogr Warning 502: Field name 'yyy xx' contains invalid characters. 'yyy_xx' will be used instead

texts:
On MSSQL tables/views that have field names with spaces or starting with a number we get the warning 502 re: replacing invalid characters with a _
As we re-injest some of these files into MSSQL and use them in other software we would like to keep the (poorly formatted) field names as otherwise we have to change every reference to the replaced fields. It could also affect AUTHORITY and Confirm etc that's used for asset management.
Is there a switch that allows us to do this?
Here's what I use
ogr2ogr -skipfailures -f "MapInfo File" -lco overwrite=yes -lco launder=no -lco relaxedFieldNameMatch=yes -lco unsetDefault=yes "D:\QG\vw_CONFIRM_PavementCondition.tab" "MSSQL:Server=xxxx;Database=GIS_Test;Trusted_Connection=True" -sql "SELECT * FROM dbo.vw_CONFIRM_PavementCondition" -a_srs EPSG:28354 -t_srs EPSG:28354 -s_srs EPSG:28354

Warning 6: dataset
D:\QG\vw_CONFIRM_PavementCondition.tab does not support layer
creation option launder
Warning 502: Field name 'Road Pavement
Maintenance Grade' contains invalid characters.
'Road_Pavement_Maintenance_Grade' will be used instead.
Warning 502:
Field name 'Road Pavement Roughness' contains invalid characters.
'Road_Pavement_Roughness' will be used instead.
Warning 502: Field
name 'Road Pavement Profile' contains invalid characters.
'Road_Pavement_Profile' will be used instead.
Warning 502: Field name
'Road Pavement Rutting' contains invalid characters.
'Road_Pavement_Rutting' will be used instead.
Warning 502: Field name
'Road Pavement Overall Condition' contains invalid characters.
'Road_Pavement_Overall_Condition' will be used instead.

Does ogr2ogr.main used within python have some additional functionality not available in the .exe? if so a solution in this would be fine as well.


texts:
How would one label at M value vertices only?

We have added M values using linear referencing and now would like to label stations at the vertices.


texts: GDAL not properly installed. Missmatch of paths or similar

texts:
I am not sure if everything is the same problem or not, but I will describe as good as I can what is going on.

In my old laptop I had QGIS installed, and when I wanted to use gdal standalone from python, I was building from source gdal2.3 and everything was working.

In my new laptop (and I freshed installed ubuntu 18.04 twice), everything has been problems. QGIS 2.18 and 3.2 don't work for a problem with libpng. (Not related I think, but something that had never happened)

Now if I try to execute some code that includes ReadAsArray in Pycharm I get the error ModuleNotFoundError: No module named '_gdal_array'.

I purged everything related with gdal (QGIS, gdal from pip, deleted the folder where I did make) and I tried to build from source again, this time doing ./configure --with-python=python3 (because i noticed that it was being installed over python2).

The most strange thing comes now: If I open a python console (which uses /usr/bin/python3.6, the same interpreter than pycharm) I can do from osgeo import gdal, but i can't import it in the script.

EDIT:

In console i can import gdal, but if I try ReadAsArray it says it can't import gdalnumeric

EDIT2:

I recompiled gdal using --prefix=/usr/lib/python3/dist-packages and added that path, but still not working...

I never had so much problems and I'm really too lost right now.


texts:
Is there any tile server that is able to render tiles for temporal data?
I need to serve tiles for rain radar data and serve them in

  tilserverurl/{x}/{y}/{z}/{time} or similar way

I've done a little research and although there are several tile servers available TileStache, mod_tile, TileCache, TileStrata,... none of them seems to support a time dimension.

Is there any open source solution available that does the job and ideally supports on demand rendering and caching?

And if there isn't any tile server that supports that out of the box, any suggestion or explanation about how to implement that or how other people implement that on production will be very helpful.

By now my best idea is to implement some homemade tileserver that is backed by N gdal datasources or mapnik maps (one for time step).
In theory is not very difficult to make it working in a simple case but make it highly performant with render queue, caching and fine tune it is another history.


texts: QGIS: Finding and tagging polygons inside a larger polygon

texts: I have two sets of polygons in shapefiles.  One set of polygons is composed of a full and exclusive set of smaller polygons.  But the smaller polygons are not tagged as which larger polygon they were used to build.  Is there a way in QGIS to determine which smaller polygons make up the larger polygons and tag the smaller polygons with the ID of the larger polygon that is belongs to?


texts:
I have a SHP point file created from CSV with QGIS 2.4.0 Chugiak. It contains some data on noise in a certain point. When I import from CSV, it seems to be fine, but when I export to SHP and open the attribute table - the real-type numbers are wrong ("reduced" number). The funny thing is, that when I open this SHP with older QGIS NextGIS 1.8 - everything is fine.

Is it some kind of bug?

I attach the printscreen:


texts: Is there a way to extract nodes that connect roads in OSM data?

texts:
Is there a way to extract only the nodes that connect links in OSM data?  I am trying to locate US Interstate exit and entrance ramp ends in order to compare with my data.  I've been able to extract end nodes but a lot of the ramps (motorway links) are a continuous line or single feature with only two ends connecting to the motorway while ignoring the secondary road.  It seems that OSM data doesn't have shared nodes but instead stacks them at connection points.  I'm going to look into extracting those but I am hoping there will be a more thorough way that will give me node attribute information if it's available.

Example is the highlighted motorway link below

EDIT:

I'm specifically looking for start and end node location of highway exit ramps (motorway_link) for all of US and Canada with the ability to differentiate between both start and end node.  Here's an example of the nodes I need.


texts: ArcGIS Server and SQL Spatial. Examples?

texts:
In-house we have an ArcGIS Server license that I am somewhat familiar with, and a license for SQL Server. I'm looking for some examples of existing web maps that use ArcGIS Server for the display of geographic data and SQL Server (Spatial) for storage. I've looked into this but haven't found any good/recent examples of this combination, especially which utilize an editing interface for modifying the stored geometries.

I'm looking for a few examples and links to documentation that will help me build a good understanding of the capabilities and requirements of these software packages before embarking on the project.


texts: DGPS Pseudorange code vs Carrier phase: What fields from the RINEX actually used?

texts:
The basic RINEX epoch measurement looks like:

Pseudorange, Carrier phase, Doppler, SNR

(In that order?)

Is this correct for post processing:

Using Pseudorange code technique (RTKLIB positioning mode DGPS/DGNSS) we can use only Pseudorange.
Using Carrier phase technique (RTK but post processing not real time) (RTKLIB positioning mode Kinematic) we can use only Carrier phase. [I assume this is not correct, since RTK using code for inital estimation and from that point improves accuracy with phase measurements]
Doppler, SNR theoretically useless for both of the above, but using them can improve results (maybe some internal algorithm taking advantage of them).

The nature of the question is to store minimum set of data for post processing.


texts: How to automate topology check/fix in QGIS?

texts:
Short version: I have a layer with 19483 features with plenty of topological errors and I would like to fix/review them as many as possible in an automated way but don't know how to do it. Additionally, since I will have to repeat the same process with many other layers with a similar amount of features and errors, I would like to create a model/script to save time.

Long version: I have this layer (7MB) with 19483 features. After using topology checker plugin I find the following topological errors:

Gaps: 228 errors. Many of these gaps are the result of polygons with very small areas that should not exist.
Must not overlap: 106 overlaps. Again, most of the problems arise because there are edges which are very close one to another (0.2mm close) but they should be in the same place instead.

Additionally, if I enable Geometry Checker plugin I found 186 sliver polygons, 159 polygons with less than 1sqm. I am not sure whether I may have even more errors that I have not yet discovered.

I am pretty new to topology (in fact, it is the first time I have to deal with it), but from what I have read, I understand that there are three ways to deal with topology checks/fixes:

Automated checks but manual fixes using Topology checker plugin.
Automated checks semi-automated fixes using Geometry Checker plugin.
Automated checks and fixes using GRASS' v.clean

If I am not wrong, after reading v.clean's documentation, it looks promising and I may solve most of the errors by using the snap tool, as if I am not wrong, would snap/move vertexes that are very close when they should be snapped instead.

Another interesting feature is that I can integrate v.clean in a model, which will make my life much better since I will have to repeat the same process with many other layers with a similar amount of features and errors.

Unfortunately, I haven't succeeded in doing so. I created the following model but I might be doing something wrong since it creates more features instead of reducing its number.

Did I do something wrong or are there better ways to deal with it?


texts:
I am using QGIS 2.18 and have a workflow that requires going offline with some of the data.  My workflow is fairly simple:

Go to Database - Offline Editing - Convert to Offline Project. Select the layers that will be edited offline.
Edit the offline layers are required.
Go to Database - Offline Editing - Synchronize

The sync process runs and is successful, but the result creates duplicate entries in the master (online) database for any of the features that were edited or created offline.

I have tried using SpatiaLite and PostGIS as my master (online, source, whatever you would like to call it) database, with the exact same results.

Any suggestions?


texts: Legend disappears when printing to PDF using QGIS?

texts:
I have a figure made in the print composer, but when I "Export as pdf, both the legend and our company logo disappear.

Any idea why?

I have not had this problem before.
Using QGIS 2.18.10


texts: Using OS data, some sea tiles are the wrong colour

texts: Using OS data and CRS ESPG:27700, Some blocks of sea tiles are showing the wrong blue colour, how do I change these tiles to the correct colour please?


texts: Which equal-area projection to use for UK?

texts:
I have a bunch of UK polygons that are saved in my PostGIS database as WGS84 projection. I would like to calculate the area of these polygons using st_area.

According to this answer I can save these polygons in any equal-area projection to do this. Is it really true that I can use any equal-area projection? For example, if I save them in US National Atlas Equal area (SRID 2163) will that work?

In any case, is there a UK-focussed equal area projection that it would be most appropriate to use?


texts: How can I reproject Shapefile in sinusoidal projection using ogr2ogr

texts:
I have a shapefile in UTM WGS84 zone 32 and I would to reproject it in sinusoidal (MODIS like).

I would try using gdal ogr2ogr command, but I am not an expert and I got no results.

I tried with:

ogr2ogr -f "ESRI Shapefile" -s_srs EPSG:32632 -t_srs "+proj=sinu" \ input.shp output.shp

with the answer:

"ERROR I: Couldn't fetch request layer output.shp"

I have also tried to insert SR-ORG:6974 instead of +proj=sinu, but this was not recognized.

May you help, please?


texts:
I'm trying to intercept SpatialLinesDataFrame (170 fields and 140000 lines) with a SpatialPolygonDataFrame using RQGIS.I'm following the instruction from github. I modified the ../processing/gui/AlgorithmExecutor.py files accordingly. My lines and polygon have numeric dataframes.

I'm running the following code:

library(maptools)
library(raster)
library(rgdal)
library(RQGIS)
setwd("/data/opera/sergio")

load("data/out/network/lines.rda") #lines
load("data/out/network/g.rda") #polygon
dir_tmp <- "/home/sergio"
my_env <- set_env()
params <- get_args_man(alg = "qgis:intersection",
                       qgis_env = my_env)
params$INPUT <- lines
params$INPUT2 <- g
params$OUTPUT <- file.path(dir_tmp, "COg.shp")

system.time(out <- run_qgis(alg = "qgis:intersection",
                params = params,
                load_output = params$OUTPUT,
                qgis_env = my_env))

Error in dirname(x) : a character vector argument expected
my_env
$root
[1] "/usr"

$qgis_prefix_path
[1] "/usr/bin/qgis"

$python_plugins
[1] "/usr/share/qgis/python/plugins"

How to fix: Error in dirname(x) : a character vector argument expected?


texts: Maplex ignores leading characters like ,,

texts:
I try to label some geologic features with Maplex. A given condition for naming is that some geological units start with a double comma ,, . So the label should look like ,,ta or ,,ya . Whenever I try to label, Maplex is not showing the leading characters. Using the standard label engine works and shows the right result, but causes more problems in conflict detection. Any ideas how to convince Maplex to give me the ,, ?

In Germany a , usually is used as decimal sign in numbers, not as separator. So my PC is set to use , as decimal sign.


texts:
I realize this isn't a geographic question as I'm not trying to make this into a geographically placed image.  But I hope you all will know how to do this.

I want to make a raster from scratch.  I have a database that I would like to be able to query and display the results as a sort of color map.  I can figure out how to do the math and raster calculating, but I can't figure out how to make a raster.  I expect that I need to add a couple of fields for the x and y coordinates, then give each record a value.  Just 1s and 0s to begin with.

I work with ArcGIS, Oasis Montaj, and Microsoft Access at work all day, but I'd like to do this in QGIS and whatever other tools are necessary.  I have spent almost no time with QGIS but would like to learn and work with it on my own projects.  But I need to point in the right direction for this one.

Thanks


texts: Does ArcGIS Server support querying multiple layers of service in single request?

texts:
I am using ArcGIS REST services from a server-side app and parsing the JSON response. I need information from several different layers on the same service, and my queries are always for points.
Is there a way to make a single request that amalgamates results from multiple layers?
At present I can meet my requirements by making multiple requests, for example to endpoints at

/arcgis/rest/services/foo/bar/MapServer/0/query?geometry=152.9259%2C-27.5516&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects...
/arcgis/rest/services/foo/bar/MapServer/2/query?...as above...
/arcgis/rest/services/foo/bar/MapServer/6/query?...as above...

My question is, am I missing some aggregated service that enables me to make a single request to multiple layers of the same service for the same spatial query?


texts: Postgis st_intersects not selecting correct features

texts:
I am running a very simple boolean selection with postgis st_intersects on two polygon tables:

select distinct on(ms.shape) ms.shape shape,ms.type_
    from riparian_corridor.merged_soils ms join riparian_corridor.riparian_open_waters o
        on st_intersects(ms.shape,o.shape);

This query returns most correct features but there are a handful of instances where the features clearly intersect and it is not selected.

The distance between the edge of that blue feature and red outline polygon is about 2 feet. Should I be using st_snaptogrid? If so, which geometry should I put it on? and should the tolerance be something around 2 feet to pick those features up?


texts:
I am trying to avoid drawing overlapping labels with vendor parameters.
I tried below option in different versions of Geoservers (2.2.1 to 2.5.2), its not working.

http://docs.geoserver.org/2.5.x/en/user/styling/sld-extensions/label-obstacles.html

My style is given below.

I have added vendor parameter in Pointsymboliser after  element.

        Point_rule1

       <PointSymbolizer>
           <Graphic>
               <ExternalGraphic>
                   <OnlineResource  xlink:type="simple" xlink:href="city.svg" />
                 <Format>image/svg+xml</Format>
                </ExternalGraphic>
              <Size>8</Size>
             </Graphic>
              <VendorOption name="labelObstacle">true</VendorOption>
         </PointSymbolizer>

 <TextSymbolizer>
          <Label>
            <ogc:PropertyName>NAME</ogc:PropertyName>
          </Label>
          <Font>
            <CssParameter name="font-family">
             <ogc:Literal>Arial</ogc:Literal>
            </CssParameter>
              <CssParameter name="fill-opacity">0.2</CssParameter>
            <CssParameter name="font-style">bold</CssParameter>
            <CssParameter name="font-size">9</CssParameter>
            <CssParameter name="font-weight">Normal</CssParameter>
          </Font>
          <LabelPlacement>
            <PointPlacement>
              <AnchorPoint>
                <AnchorPointX><ogc:Literal>0.0</ogc:Literal></AnchorPointX>
                <AnchorPointY><ogc:Literal>0.0</ogc:Literal></AnchorPointY>
              </AnchorPoint>
             <Displacement>
                <DisplacementX><ogc:Literal>8</ogc:Literal></DisplacementX>
                <DisplacementY><ogc:Literal>0</ogc:Literal></DisplacementY>
              </Displacement>
              <Rotation><ogc:Literal>0</ogc:Literal></Rotation>
            </PointPlacement>
          </LabelPlacement>
          <Halo>
            <Radius><ogc:Literal>1</ogc:Literal></Radius>
            <Fill>
              <CssParameter name="fill">#FFFFFF</CssParameter>
              <CssParameter name="fill-opacity">1</CssParameter>
            </Fill>
          </Halo>
          <Fill><CssParameter name="fill">#333333</CssParameter></Fill>

         <VendorOption name="goodnessOfFit">1.2</VendorOption>
         <VendorOption name="removeOverlaps">true</VendorOption>
         <VendorOption name="labelObstacle">true</VendorOption>
        </TextSymbolizer>

      </Rule>
    </FeatureTypeStyle>
</UserStyle>

Anyone implemented this vendor option successfully?


texts: Get rid of overhang beyond masking border in ggmap?

texts:
I built a terrain map of Nevada and plot the state border around it. The resulting map has a buffer-like overhang across the state border which I used to mask the raster. How can I get rid of the overhanging data? The code below makes this example easily reproducible.
    library(ggmap)
    library(raster)
    library(sf)
    library(ggplot)
    library(USAboundaries)
    state <- us_states(states = c("Nevada"))
    map3 <- get_googlemap(center = c(lon=-116.888369, lat=39.509636), zoom = 4, maptype = "terrain", color= "color",
                          style = "feature:road|visibility:off&style=feature:administrative|element:labels|visibility:off&style=feature:administrative|element:geometry|visibility:off&style=feature:landscape|visibility:on|saturation:-100|lightness:-10&style=feature:landscape|element:labels|visibility:off&style=feature:poi|visibility:off&style=feature:transit|visibility:off")
    map.rast <- ggmap_rast(map = map3)
    map.only <- mask(map.rast, state)
    basemap.df <- data.frame(rasterToPoints(map.only))

    ggplot() +
      geom_point(data= basemap.df, aes(x=x, y=y,col=rgb(layer.1/255, layer.2/255, layer.3/255))) +
      scale_color_identity() +
      geom_sf(data=state, color= "black", fill=NA, size=0.3, inherit.aes = F) +
      coord_sf()


texts: Exporting pandas dataframe as shapefile using ArcPy gives ERROR 000732

texts:
I'm trying to export my dataframe as a CSV file and as a shapefile. I'm able to export it as a CSV but get the following error on the last line:
ERROR 000732: Input Features: Dataset in_memory\x37C54D46_4BEE_4D63_B0B5_7EF0B54FEE74x does not exist or is not supported
Failed to execute (FeatureClassToFeatureClass).

Any idea why it doesn't recognize my dataframe?
Is there another way I can export it as a shapefile?
outDir = r"C:\Users\Test"
now = datetime.datetime.now()
currentDate = str(now.month) + "_" + str(now.day) + "_" + str(now.year)
sdf = pd.DataFrame.spatial.from_featureclass(fc)

#Add fields to point
sdf['NADatum'] = 'WGS84'
sdf['Shape Type'] = 'Point'
sdf['Method'] = 'GPS'

#Pull only the necessary fields
sdf = sdf[['REPORTER','RECOWNER','SCINAME','COMNAME','RECBASIS','GlobalID','OCCSTATUS','MANAGESTAT','POPSTAT','OBSDATE','COUNTRY','STATE','COUNTY','PRIVATE','ESRIGNSS_LATITUDE','ESRIGNSS_LONGITUDE','NADatum','Method','Shape Type','ABUNDANCE','INFESTAREA','INFESTUNIT','PCTCOVER','VISITTYPE','COMMENTS']]

sdf.columns =['Reporter','ReportOwnershipID','ScientificName','CommonName','Recordbasis','UUID','OccurenceStatus','ManagementStatus','PopulationStatus','ObservationDate','Country','State','County','Private','Latitude_Decimal','Longitude_Decimal','NADatum','Method','ShapeType','Abundance','InfestedArea','InfestedAreaUnits','PercentCover','VisitType','Comments']
#Export point as CSV
outputCSV = os.path.join(outDir, "Points_" + currentDate + ".csv")
sdf.to_csv(outputCSV, index = False, header=True)
outFeatureClass = "Points_"
arcpy.FeatureClassToFeatureClass_conversion(sdf, outDir,outFeatureClass)


texts: I've created two different buffers (1 mile and 3 miles) in ArcMap, I need to see how many people live within 1 and 3 miles. In other words, I need to subtract the bigger buffer from the smaller one. I have no idea how to do it.


texts: Summarising all fields in attribute table using ArcPy?

texts:
I was hoping to run a script that will iterate through all fields in an attribute table, SUM each field and produce the result in a table.

The script below will essentially do what I want with one field but is there a way for it to iterate through all fields?

import arcpy

from arcpy import env

env.workspace = "C:/data/Habitat_Analysis.gdb"

arcpy.Statistics_analysis("FIELD_NAME", "C:/output/output.gdb/TABLE_OUTPUT", [["SUM"]])

(the above sample was sourced here: http://resources.arcgis.com/en/help/main/10.1/index.html#//00080000001z000000)


texts: How can I transfer the values of a grid of smaller polygons to another grid with larger polygons in ArcGIS?

texts:
I have a population density grid containing 100x100m polygons. I would like to transfer the density values into an existing overlaying grid of larger polygons (250x250m). This might of course include modifications of the actual values.

any ideas how to do that?


texts: How to georeference to the US base in Guantanamo?

texts: I make my first steps on CartoDB. I would like to assign data to the US enclave at Guantanamo in Cuba and I can't. Yet the CartoDB base map shows the enclave.


texts: Compare two fields not equal to each other in a join from a feature layer

texts:
I am wanting to find features not equal to each other in a feature layer in Python, and I am not finding the right SQL to use. The two features are joined to each other on a MasterID field. When using Select By Attributes query TLNO_GIS <> TLNO, the results are successful. When running the same query in Python (below), the error message returns:

  ERROR 000358: Invalid expression ("TLNO_GIS" <> "TLNO")
  Failed to execute (MakeFeatureLayer).

arcpy.MakeFeatureLayer_management(Production_ISGIS_COH_Address, ISGIS_COH_Address_LayerTaxlot, "(\"TLNO_GIS\" <> \"TLNO\")", "", "OBJECTID OBJECTID VISIBLE NONE;MasterID MasterID VISIBLE NONE;TLNO_GIS TLNO_GIS VISIBLE NONE;TLNO_Taxmap TLNO_Taxmap VISIBLE NONE;TLNO_status TLNO_status VISIBLE NONE;SiteAddress SiteAddress VISIBLE NONE;SitusNum SitusNum VISIBLE NONE;StreetPrefix StreetPrefix VISIBLE NONE;StreetName StreetName VISIBLE NONE;StreetSuffix StreetSuffix VISIBLE NONE;Unit Unit VISIBLE NONE;UnitNum UnitNum VISIBLE NONE;Building Building VISIBLE NONE;BuildingNum BuildingNum VISIBLE NONE;PO_City PO_City VISIBLE NONE;PO_State PO_State VISIBLE NONE;PO_Zip PO_Zip VISIBLE NONE;AddrType AddrType VISIBLE NONE;Description Description VISIBLE NONE;Status Status VISIBLE NONE;AddrJur AddrJur VISIBLE NONE;Office Office VISIBLE NONE;Floor Floor VISIBLE NONE;BuildingID BuildingID VISIBLE NONE;RegionID RegionID VISIBLE NONE;CreateID CreateID VISIBLE NONE;CreateDate CreateDate VISIBLE NONE;Editor Editor VISIBLE NONE;EditDate EditDate VISIBLE NONE;SubPudComplexName SubPudComplexName VISIBLE NONE;PlatLotNum PlatLotNum VISIBLE NONE;PlatRecDate PlatRecDate VISIBLE NONE;X_Coord X_Coord VISIBLE NONE;Y_Coord Y_Coord VISIBLE NONE;Angle Angle VISIBLE NONE;Shape Shape VISIBLE NONE;GlobalID GlobalID VISIBLE NONE;FromSitusNum FromSitusNum VISIBLE NONE;Phase Phase VISIBLE NONE;MasterID MasterID VISIBLE NONE;TLNO TLNO VISIBLE NONE")


texts: Using ArcGIS API for JavaScript 3,X Directions widget?

texts:
Can I select only one country in the Directions widget in ArcGIS API for JavaScript 3,X?

For example: I want to select only UK. If the user writes Spain he will find nothing.


texts: Getting ST_IsSimple reason and detail similar to ST_IsValidDetail in PostGIS?

texts:
How can you get ST_IsSimple reason and detail similar to ST_IsValidDetail in PostGIS?
When validating geometries I can use these handy functions to not only check if geometry is valid, but also the detail when not:
SELECT reason(ST_IsValidDetail(geom)) AS reason, ST_AsText(location(ST_IsValidDetail(geom))) AS location FROM my_polylines WHERE ST_IsValid(geom) = false;

This returns reason and location of the problem. I would like to get similar information for geometries where ST_IsSimple(geom) = false. Are there any PostGIS functions to achieve that? Or any guide how to achieve it through combination of PostGIS functions?
From "Ensuring OpenGIS compliancy of geometries":
A LINESTRING is simple if it does not pass through the same POINT twice (except for the endpoints, in which case it is referred to as a linear ring and additionally considered closed).
So talking about LINESTRING only, there is always only one reason - 'self-intersection' (including self-overlap). The goal is to find out how many and where exactly.


texts: Openlayers HiDPI/Retina Screen with standard resolution raster tiles?

texts:
My Openlayers application can display tiled rasters from arbitrary online sources.  It will parse capabilities from any user-specified URL for WMTS, TMS, OSM, XYZ or WMS services and present the user with a list of available layers which they can choose from.

Many of these layers are only available in standard resolution, but most of the users will be using Retina or high DPI screens.  The result is that the map features are all too small, thin and spidery looking and text is very difficult to read for people with standard eye sight (impossible for anyone with aging eyes).

Is there some way to get Openlayers to display these standard resolution tiles at double the resolution (or x1.3 or x3 or whatever the screen scale is)?

Ie, To display the normal 256x256 tiles at 512x512 pixels?

I know this would make the rasters look a little blurry, but in theory no more so than if they were displayed on a standard resolution screen.  At least they'd be big enough to view/read properly - and I can give the user the option of whether to use this feature or not.  I've verified that the rasters look fine at double size by magnifying the view in which they are displayed by 2x, but of course this magnifies everything else in the web page too, which is undesirable  Eg, widgets are double-sized, vector layer symbology gets blurry.

I can easily scale vector symbology, just my multiplying the various width and size parameters by the screen scale factor (this is tested and working), but I cannot do anything similar for rasters.  So I think the best work around would be to show them at double size.  Eg, grab the tiles that would normally be used if zoomed out by 2x and display them at 2x their actual size.

It seems that this should be possible but I can't figure out how it would be done.  Perhaps redefining all the layer's resolutions and setting tilePixelRatio = 0.5 might be part of the solution, but my Openlayers knowledge is limited and I'm not sure if this is even on the right track.

I expected a lot of other people must be trying to deal with the same issue, but I've been unable to find any information that fits this scenario.

Is this possible?  How can it be done?

JSFIDDLE UNSUITABLE TO DEMO THIS ISSUE:

I have created a jsfiddle to demonstrate this ( https://jsfiddle.net/nsands/n3cbwr95/20/ ) but jsfiddle is not suitable in this case, as it and/or iOS does its own strange scaling on the iPhone retina screen.  With iPhone in portrait mode, everything is too tiny (not just the map, but even the HTML and Javascript panes have text that is unusually small), and when in landscape mode, the map text is large enough, but appears to have been scaled up (it is blurry - actually the kind of behaviour I WANT to achieve elsewhere, but only for raster layers).

ACTUAL DEMO OF THIS ISSUE:

However, this issue can be demonstrated using the WMTS example on the Openlayers examples page at:  https://openlayers.org/en/latest/examples/wmts.html

I've loaded this up on both a standard resolution desktop PC and on a 'retina' (HiDPI) iPhone 6S and focused on the same location at the same zoom level on both devices.  Below is what it looks like on both (note that the retina image is actually double the resolution shown here, but is the same size as shown here - and it is the size of the elements that is the problem for me).

Note that the elements in the HiDPI/'Retina' version are about half the size of the elements in the standard resolution version.  For some people, this is too small to read, and on some maps the elements are rendered even smaller than this.

OpenStreetMap uses fairly large fonts when rendering its tiles.  Some map services (such as that provided by my state government) use smaller fonts, and are completely unreadable in Openlayers on retina/HiDPI screens (ie, https://services.thelist.tas.gov.au/arcgis/rest/services/Basemaps/Topographic/ImageServer/WMTS/1.0.0/WMTSCapabilities.xml ).

STANDARD RESOLUTION:

'RETINA' (SCALED TO ACTUAL SIZE):


texts: Can an arcpy.da.SearchCursor be used on a dataset residing in a ParcelFabric?

texts:
I'm using ESRI's LocalGovernmentMI.gdb data model and am managaing parcel data in a ParcelFabric dataset.  I'm trying to use a search cursor on the ParcelFabric_Parcels dataset but I recieve this ...

  RuntimeError: cannot open 'P:\Mapping\Data\LocalGovernmentMI.gdb\ParcelEditing\ParcelFabric\ParcelFabric_Parcels'

I'm using this but of code ...

import arcpy

fc = r'P:\Mapping\Data\LocalGovernmentMI.gdb\ParcelEditing\ParcelFabric\ParcelFabric_Parcels'
field = "Name"

with arcpy.da.SearchCursor(fc, field) as cursor:
    for row in cursor:
        print (row)

If I export the dataset and run the search cursor on that data it works.  Any thoughts?


texts: Adjust polygons by edition (freehand and reshape tools used together)

texts:
I would like to adjust the polygons of my vector file. This vector is a land use classification with classes of vegetation, human uses, etc. So the form of the features classes are so irregular and classification isn't totally correct. So, working only with reshape is so hard to digitalized adjusting, many clicks in order to draw the form feature of image layer perfectly.

Is there possibility to use reshape and freehand together qgis?

Is there any suggestion to make this kind of edition more easy and fast?


texts: Export all data driven pages to PDF using Python script

texts:
I am a novice Python user but I am trying to create a script which will allow me to export all of the data driven pages within an mxd to a single pdf file (without having to open the mxd). I currently have the script below which works for exporting single pages. And the highlighted section was my initial (poor) attempt to export all pages based on what I have seen or read.

# Python script to run in windows.
# 1) Copy script into folder to run script
# 2) Edit quality and dpi variables if desired

import arcpy, os

# get current directory python script is in
cwd = os.getcwd()

# quality to export, BEST, BETTER, NORMAL, FASTER, FASTEST
quality = "BEST"
# dots per inch resolution
dpi = "300"
res = int(dpi)

# Ask user to confirm settings and directory
print "Export ArcMap MXD(s) to PDF(s), " + quality + " Quality, " + dpi + "dpi to " + cwd + "?"

# pauses for user to press key to continue
os.system('pause')

# set folder path to current directory
folderPath = cwd

# prints message to console for user
print "IN PROGRESS...PLEASE WAIT..."

# counter to count number of PDFs
n = 0

# create fullpath from filename and folder path
for filename in os.listdir(folderPath):
    fullpath = os.path.join(folderPath, filename)
# extract file extension
    if os.path.isfile(fullpath):
        basename, extension = os.path.splitext(fullpath)
        if extension.lower() == ".mxd":
# advance counter to count PDFs
            n = n + 1
            mxd = arcpy.mapping.MapDocument(fullpath)
# export tool with parameters
            arcpy.mapping.ExportToPDF(mxd, basename + '.pdf', resolution=res, image_quality=quality, colorspace="RGB", picture_symbol="VECTORIZE_BITMAP", layers_attributes="NONE")
            print basename + "...exported"

            del mxd

# prints message to console for user
nstr = str(n)
print nstr + " pdf(s) exported to " + folderPath
print "DONE"

# pauses for user to press key to continue
os.system('pause')


texts: Finding area of individual multipart parts using PyQGIS

texts:
Im trying to find the area of each part in multipart polygons. I know some polygons have multiple parts, so why cant I access the individual polygon parts like this (why is parts list empty):
lyr = iface.activeLayer()
f = next(lyr.getSelectedFeatures())
parts = [p for p in f.geometry().parts()]

print(len(f.geometry().asMultiPolygon()))
print(parts)

2
[]


texts: Why does textAlign in Openlayers look opposite?

texts:
As far as I know with CSS knowledge, textAlign:left means the text should be stick left of the element. W3School textalign

but in Openlayers, it looks working opposite.
Openlayers vector label example.
click the link above, and choose left from align select box, then press refresh then the texts moves to right side like textAlign:right.

How does textAlign work in OpenLayers?


texts: Creating a weighted rating model/script for vector data

texts:
I'm using ArcMap and trying to create a model or script that can rate routes based on a weighted calculation of attributes. Right now the individual attributes don't have rating values which I could weight and calculate to a total rating value.

I have a feature that has several attributes. Each attribute has a set of domains that define possible values in each field.

I also have a separate table that has data and is related to each feature through a relationship class (route obstructions). The values in this table area also defined from a set of domains.

What I would like to do is using model builder (which I prefer) or a script,
calculate a numeric value into a new field to form a weighted, calculated value, based on the attributes of the feature.

All these fields have different domains, and essentially what I want is to assign each option from the domains a specific value so that a new field can be calculated based on these values.

The feature also has a related table, and I would like to use the values in the fields of this table in the calculation as well.


texts:
I have a table with locations and labeled data and I need to apply Delaunay Triangulation on specific labeled rows.

For example, I have some locations labeled as schools, hostels, gardens, ... and I need to apply triangulation on those specifically labeled items separately.

To solve this I can create as many new tables as labels of interest, but even if I could create such tables it would be daunting task and clearly not very intelligent approach.

Perhaps I could write Python script, but first I thought to ask if there is simpler solution, either in ArcMap or QGIS.


texts: Incomplete Data - Missing line segments from Oracle data

texts:
I have setup GeoServer 2.5.2 on Windows.  Everything was going smoothly until I actually had to display my geometry data from an Oracle Database Table.  Some lines are rendered but some are not rendered.  The issue seems to be based on the type of lines (curves, etc). This issue does not happen if exporting from oracle to SqLite or Global Mapper.  In that case, all lines are kept as in the original data from Oracle.  The issue seems to be most likely with the oracle plugin (geoserver-2.5.2-oracle-plugin.zip).

Has anyone else ran into this issue and if so, how did you manage to render everything coming from Oracle in GeoServer?


texts: Finding the date of acquisition of Google maps in QGIS

texts: How do i find the date of acquisition of the Google map on QGIS? Please help.


texts: SSDM : Error in rasters stack and SSDM's model prediction

texts:
I'm actually working on Environemental data to build an SDM model.
I faced some problems during data processing.
 I extracted my environnemtal  from ENVIREM database, res=30 arc-sec.
 I checked my rasters before stacking them, there was 0 NA's.
When I stack them I get a lot of NA values as shown below:

AnPETM_1  topwetM_1 ArdindexM_1  PETSeasM_1 PETWettQuartM_1
Min.    2.669464e+01     885.400       6.17       63.69    1682.990           25.81
1st Qu. 4.269421e+01    1333.732      10.54       78.82    3789.467           61.24
Median  5.209351e+01    1464.990      12.32       92.96    4974.755           88.65
3rd Qu. 6.036992e+01    1640.975      13.51       97.47    5840.647          115.53
Max.    8.304529e+01    2012.200      15.17       98.90    7386.350          213.46
NA's    2.594705e+06 2594248.000 2594248.00  2594248.00 2594248.000

Is there any solution to keep real values in my rasterstack during prediction process.

2/ I used my raster stack in SSDM for single species distribution.
I got this;

 ***Warning messages:
1: In min(pred[which(obs == 1)]) :
  no non-missing arguments to min; returning Inf
2: In min(pred[which(obs == 1)]) :
  no non-missing arguments to min; returning Inf
3: In cor(predicted.values, o.predicted.values) : l'écart type est nulle
4: In cor(predicted.values, o.predicted.values) : l'écart type est nulle
5: In cor(predicted.values, o.predicted.values) : l'écart type est nulle
6: In cor(predicted.values, o.predicted.values) : l'écart type est nulle
7: In cor(predicted.values, o.predicted.values) : l'écart type est nulle***


texts: Parcel Owner to see only his parcel in map application

texts:
How can i create a map where user can only see his parcel. parcel are defined has one unique ID. I have Arcgis online Server and Desktop Advanced license.
the map to be published in a web application which will show owners details reading from SQL server database.


texts:
Is it possible to use a variable in a layer filter? For example, I'd like to have

"STATE" = @current_state

to only display the features within a certain state. Trying to do this gives me an SQL syntax error.


texts: Leaflet with GeoServer in background too slow for GeoTiff

texts:
I am implementing a Leaflet web-map application to display several GeoTiffs. The GeoTiffs are each ~300MB and published by GeoServer as WMS. In Leaflet I add them through L.tileLayer.wms.

I also perform some grouping of the layer before adding them to the map.

My problem is now that the layers published by the GeoServer load very slow in the Leaflet application, whereas the OpenLayers preview offered by GeoServer is pretty fast.

Any Idea what this might cause?


texts: Is there a more direct approach to time series correlation in the Google Earth Engine?

texts: I want to know how the time attributes of ee.Image and ee.ImageCollectionobjects exist in Google Earth Engine. I am looking for whether there are operations on time series in GEE, such as resampling, sliding window, etc. The bigger thing is to know if the Python Api can be combined with other libraries, such as pandas.


texts:
I want to export a lot of raster layers as rendered images. After looking around, I got to this code:
myDir = 'e:/Data/'
layers = [layer for layer in QgsProject.instance().mapLayers().values()]
pipe = QgsRasterPipe()
for layer in layers:
   extent = layer.extent()
   width, height = layer.width(), layer.height()
   renderer = layer.renderer()
   provider=layer.dataProvider()
   crs = layer.crs().toWkt()
   pipe.set(provider.clone())
   pipe.set(renderer.clone())
   file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
   file_writer.writeRaster(pipe,
                           width,
                           height,
                           extent,
                           layer.crs())

But when I run it, I get this error:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 6, in <module>
AttributeError: 'QgsVectorLayer' object has no attribute 'width'

I'm using QGIS 3.12.2.
Any suggestions?


texts: Central deployment of startup.py in QGIS

texts:
In QGIS 3 on Win10 we use a startup.py script for setting up e.g. database connections and a couple of other things, residing in C:\Users\me\AppData\Roaming\QGIS\QGIS3

What we are after is to centrally deploy the startup.py without having to distribute it to each of our clients at every script change.

First try was to make use of the PYQGIS_STARTUP environment variable, but this fails, because startup.py makes use of qgis.* modules which are of course not available before QGIS initialization is complete.

How can we centrally deploy our startup.py in order to apply changes on every client without having to distribute it?


texts: Georeferencing a raster with known Lon/Lat center coordinate

texts:
I would like to digitize a scanned map, see example below.

The projection and datum information is available on the map. I am familiar with georeferencing of a scanned image that contains multiple latitude/longitude coordinates.

Unfortunately, this particular map has only one lat/long coordinate, map's geographic centre, plus distance measurement (i.e., in yards). I am thinking maybe I could convert the distance measurement into degrees, but I don't know how. Is this possible?


texts: Using Aggregate and Infowindow in CartoDB

texts:
I'd like to use CartoDb to map historic  addresses of a city (link to demo). Individual addresses are mapped as points on the map. Detailed information for each address point should be displayed in the infowindow The geocoded addresses are stored like this:

NAME   PROFESSION ADDRESS ... (+ CartoDb fields like cartodb_id, the_geom_webmercator)
Smith  shopkeeker grand-rue 10
Martin baker      grand-rue 11
Wagner grocery    grand-rue 11

As long as there is only one person living at an address, I get the desired infowindow by simply doing a select * from _table_name_. When two or more persons are associated with the same address I obviously would need some aggregate function in PostgreSQL.

I am experimenting with:

SELECT  address, string_agg(name || ', ' || profession, '; ' ORDER BY name, profession) FROM table_name group by address

I am running into two problems:

In order to map, CartoDb expects the columns 'cartodb_id' and 'the_geo_webmercator' to be included in the SELECT statement. When doing so, the result of the query is not grouped by 'address' in the desired way; instead the unique 'cartodb_id' is used for grouping. So, the following is not grouping according to 'address':

SELECT  address, string_agg(name || ', ' || profession, '; ' ORDER BY name, profession), the_geo_webmercator, cartodb_id FROM table_name group by address, the_geo_webmercator, cartodb_id

The output of 'string_agg' is treated as one field in the design for the infowindow.

Grand-rue 11
Martin, baker; Wagner, grocery

What I would like to have, though, is a separate treatment of 'name' and 'profession' in the infowindow, like:

Grand-rue 11
1. Martin
   baker
2. Wagner
   grocery


texts: Where can I find the highest resolution shapefiles for counties or states

texts: I am looking for a download source which has the highest resolution (greatest detail/accuracy) for all counties in the US. More specifically I run into this problem mainly when dealing with the shorelines. I realize that one can download high resolution shoreline data, overlay the tiger shapefiles for the counties and then clip them but I figured someone might have already done this.


texts: Converting Service Area polygons to raster using ArcGIS Network Analyst?

texts:
I have just made some very nice service area polygons and would like to convert these to rasters for further analysis. When I try polygon to raster (field: to break), it results in an empty raster. I thought maybe there is a problem because there are so many different polygons... so I attempted to dissolve the polygons according to break values, but the resulting polygon shapes are totally different from the ones in the service area analysis. In other words, the shapes from the network analysis are not preserved.

Why?

Update: Here is a screenshot of the map and assoc. attribute table.

TO respond to one comment: I have tried to export them as polygons, but when converting these to raster, the shape is not preserved, assuming because there are so many polygons involved. Wish there was a simple fix to this for quick analysis of network polygons and service area findings. I have uploaded a screenshot the raster result as well.


texts: Zonal statistics - computation time issue with assigning zones

texts:
I would like to upscale a 1/360 degree resolution raster (64800 x 129600) to a 1/12 resolution raster (2160 x 4320) by assigning to each cell in the low-res raster the number of cells in the high-res raster that are contained in that cell AND that have a certain value.
In order to use the 'zonal' function, I first need to assign zones to the grid points of the high-res raster. Every 30x30 block of grid points should be one zone. (I.e. the zone raster should equal '1' in the top left 30x30 grid points, and 2160*4320=9331200 in the bottom right 30x30 grid points).
Here is the code I am using to generate this zone raster ('raster_high' is the original high-res raster.)

zones <- raster_high;
for(i in 1:64800)
{
  j <- floor(i/30);
  zones[i,] <- rep(4320*j+1:4320*(j+1), each=30);
}

However, even only the first step i=1 takes ages. Is there a quicker way to do this?


texts: Locate maximum cell within a radius r around a central cell

texts: I have a raster. Is it possible to find the location of a maximum cell within a defined radius'r' around a central cell.


texts: How to determine number of features in GeoJSON file using Python or JavaScript?

texts:
Is there a best way to determine the number of features in a GeoJSON file?

Ideally using Python or JavaScript.


texts: Converting google map api coordinates into meters

texts:
How would I go about converting the results of my map bounding box into meters within javascript?

For example, if I have the coordinates 51.52366625446907, -0.07467269897460938 decimal degrees, can I project that into something like 555000 155000 meters?

(If it is easier, I could try to export my raster from ArcGIS in decimal degree format; I don't know which option would be better).


texts: Overpass query: Search range of integers not working for 2 digit numbers

texts:
I would like to search for all admin levels with the values between 7 to 11 using the following query (http://overpass-turbo.eu/s/Dn5):

(relation[admin_level~"^[7-11]$"]["wikidata"]({{bbox}});

The error says:

  static error: Invalid regular expression: "^[7-11]$"

Surprisingly a search query with ^[3-6]$ works perfectly fine. So it seems to be related to the 2 digit number.

Does anybody know how to improve this query?


texts: What is the easiest way to create map using HDF-5 data subset?

texts:
I need to visualize some area using values provided in HDF-EOS version 5 file and create map using those values and other GIS layers.
Unfortunately, my QGIS do not recognize the HDF5 file.

Do I really need to compile GDAL manually as described at: SciGeo?

Maybe there's an easier way to do that? I thought about conversion to the GeoTIFF. However HEG tool crashes during conversion.

I opened my product in HDFView, but I don't see any option to generate GeoTIFF in this application.

Last, but not least, I've tried R!
I read appropriate data set, but I don't know what to do next.
Is there any possibility to save this data in different format?
Or maybe I am able to (how? I can't find this in h5 package documentation)

My current code is:

pacman::p_load(h5)

# download HDF-5 file
destfile <- 'SMAP.h5'
if(!file.exists(destfile))
{
  download.file('ftp://n5eil01u.ecs.nsidc.org/SAN/SMAP/SPL2SMAP.003/2015.07.07/SMAP_L2_SM_AP_02293_D_20150707T084111_R13080_001.h5',destfile)
}

# read file
file <- h5file(destfile ,'r')
dataset <- file["Soil_Moisture_Retrieval_Data"]["soil_moisture"]

# @Spacedman: Does your R have a GDAL that supports these files?
subset(rgdal::gdalDrivers(),grepl("HDF",long_name))
# it outputs:
#         name    long_name create  copy
# 51 HDF4Image HDF4 Dataset   TRUE FALSE
# 53 HDF5Image HDF5 Dataset  FALSE FALSE


texts: How to filter OSM files based on tags using Osmosis?

texts:
I am working with .osm files and am trying to filter entities based on tags. I need to view only entities containing both source and attribution tags. I used the following to successfully filter all source tags only:

osmosis --read-xml file = "northern_ireland.osm" --way-key-value keyValueList="source"
--used node --write-xml northern_ireland_source.osm

From another tutorial, I've noticed that I would need to use a pipe: http://wiki.openstreetmap.org/wiki/Osmosis/Detailed_Usage_0.38

So I did the following:

osmosis --rx northern_ireland.osm --accept-ways source=* outPipe.0 = SOURCE

However, upon trying such code, I get a syntax error stating:

Argument does not contain a name before the '=' (i.e. name = value).

Any suggestions?