texts: Get signal on delete group on QGIS 2.18.16

texts:
I tried :

QgsProject.instance().layerTreeRoot().willRemoveChildren.connect(self.onWillRemoveChildren)

It works if I delete a group within in the root node. But if I delete a layer within a group QGIS crashes (core dump).

my code :

class Project:
    def __init__(self):
        self.name = ""
        self.layer_tree_group = self.set_layer_tree_group()
        self.layers = []

    def set_layer_tree_group(self):
        self.layer_tree_group = QgsLayerTreeGroup()
        self.layer_tree_group.setName(self.name)
        self.layer_tree_group.setCustomProperty("Project", self.project.name )

    def add_layer(self):
        point_layer = QgsVectorLayer("Point?crs=4326", "project_layer", "memory")
        QgsMapLayerRegistry.instance().addMapLayer(point_layer, False)
        self.layer_tree_group.insertLayer(0, layer)

 class PluginInterface:

    list_project = []

    def __init__(self):
        self.active_project = self.group_selected
        QgsProject.instance().layerTreeRoot().willRemoveChildren.connect(self.onWillRemoveChildren)

    def onWillRemoveChildren(self, node, indexFrom, indexTo):
        print "Will remove", node.name(), indexFrom, indexTo

The aim is to delete instance of Project if a user remove a group or a layer in the legend.

How can achieve it ?


texts: How to style the WFS layer in openscales?

texts:
Possible Duplicate:
  Styling the WFS point layer (openscales 1.1 & Flex 3)

I am very new to openscales...
I want to change the styles of the WFS point layer dynamically by giving rules. I want to show some jpeg/image as the point in the openscales.

Any idea how to define the rules to change the styles of default point style.


texts: Why are the layer-tree-layer source attributes absolute even if my project is set to relative paths?

texts:
When I save a group as a layer definition file (both from the context menu and via the python api) the
layer-tree-layer source attribute is incorrectly absolute:
    source="/Volumes/gis/master/Raster/filename.tif"

while the maplayer datasource is correctly relative:
    <datasource>./Raster/filename.tif</datasource>

If I change the .qlr text file (of the layer definition) to set the layer-tree-layer source attribute to relative, the layer definition will fail to load.

If I change the .qgs xml text file (the project file) to set the layer-tree-layer source attribute to relative, the group will show up empty  (failing silently) when loading the project file.

Is this a bug?


texts: Trouble installing Rasterio module for python 2.7 on windows

texts:
Can anyone please explain to me in step by step instructions how to install rasterio with Python 2.7 on Windows. I have Windows 8 and have tried information on Link

I have installed all the dependencies mentioned and when I try to run the following command on command prompt:

$ rasterio-0.24.1-cp27-none-win_amd64.whl

It displays an error:

rasterio-0.24.1-cp27-none-win_amd64 is not a supported wheel on this platform.

I have tried all the binaries given at http://www.lfd.uci.edu/~gohlke/pythonlibs/#rasterio but none has worked.

It's installation requires GDAL also and when I try any of the GDAL binaries it gives the same error message as that for rasterio but nevertheless I have installed GDAL with  OSGeo4W installer. Now how do I go about installing rasterio? Please suggest something.


texts: Creating gradient point on top of shapefile using QGIS?

texts:
I am attempting to create a point that scales in relative size to a particular integer (frequency) and is gradiented to a particular integer value, and then overlay that point onto each specific polygon of a shapefile (USA)

The point should look something like this where darker colored points denote higher numbers and larger circles denote higher frequencies.

I have frequencies all values (included frequencies) joined to the USA shapefile, but am unable to create the points.

Are they buffers?

I'm new at this (QGIS).


texts: Create a raster dataset containing the number of returns from a las file?

texts:
I already created DEMs, DTMs, DSMs from a LiDAR (LAS) file with a little help from the internet. However, right now I'm kind of stuck in my research.

I want to create a raster file which contains for each pixel the number of returns from the las file. So that the result would be a raster file with 5 classes (values) because the LiDAR system used to generate the data was able to differentiate in up to 5 LiDAR returns from one pulse.

ArcGIS 10.2.2, ERDAS Imagine 2014 and LAStools are available to do that job.
I would be very thankful for any kind of help!


texts:
I use the pitfree algorithm to create the CHM of 89 .LAS files.

chm = grid_canopy(lasnorm, 0.5, pitfree(thresholds = c(0,2,5,10,15,20), subcircle = 0.2))

In the entire code, I go trough each files to produce a segmentation raster with some metrics. Depending on the settings I use for the pitfree algorithm, I will get an error on different LAS files.

I tried to change the parameters (resolution, subcircle, thresholds), but each time I would end out having that same error ("Not a matrix"), although, not for the same .LAS. I did not play much with the max_edge argument because I am not yet entirely sure of how it affects the algorithm.

The error I get from running the code in a loop (tSearch()) seems to refer to the triangulation part of the code.

Here's the error I would get running the code with a loop

Error in tSearch(D, P, X, threads) : Not a matrix.

Here's the error I would get running the code with LAScatalog

 An error occurred when processing the chunk 7. Try to load this chunk with:
 chunk <- readRDS("...\AppData\Local\Temp\RtmpQ3rscv/chunk7.rds")
 las <- readLAS(chunk)
 Not a matrix.

I would like to know what I should look into to avoid this error but keep the same parameters for all my .LAS?

As suggested, I tried manipulating the chunk from the error. I put here my results in case it can be useful.

I am able to run pitfree() on it with different parameters (las is the chunk):

chm18 = grid_canopy(las, 0.5, pitfree(c(0,2,5,10,15,20))) #Without subcircle (error)
Error in tSearch(D, P, X, threads) : Not a matrix.

chm18 = grid_canopy(las, 0.5, pitfree(c(0,2,5,10,15,20), subcircle = 0.2)) #with subcircle (works)
chm18 = grid_canopy(las, 0.25, pitfree(c(0,2,5,10,15,20))) #res to 0.25 (works)
chm18 = grid_canopy(las, 1, pitfree(c(0,2,5,10,15,20))) #res to 1 (works)

I tried with the edge parameters too.

chm18 = grid_canopy(las, 0.5, pitfree(c(0,2,5,10,15,20), c(0,1)))
Error in tSearch(D, P, X, threads) : Not a matrix.

c(0,3), c(2,3), c(0,4) all worked though

Lastly, I changed the thresholds values :

chm19 = grid_canopy(las3, 0.5, pitfree(c(0,2,5,10))) #removed the 2 last tresholds (works)

chm19 = grid_canopy(las3, 0.5, pitfree(c(0,2,5,10,14))) #tried last threshold set other than 15 (works)

Any threshold other that 15 would work (10,11,12,13,14,16,17,18,19,20)


texts: Is there an equal area basemap for the Alaska region of North America?

texts:
I am looking to apply an equal area basemap to my data. Although I do not receive an error when applying the data in the same arcmap document, the data is being spaced in an undesirable fashion. The spacing of the point data is not wrong, I simply want each pixel to represent the same amount of area on the map. To be clear, the points are indeed representing the same area but it is the visual representation of the area which is being affected.

Data without basemap, hence no distortion/seperation between pixels:

Data with basemap and distortion between pixels:


texts:
I use PostgreSQL 9.2 but geometry_column table not automatically exist like spatial_ref_sys table, but if I use PostgreSQL 9.1 it's exist.
How to make geometry_column automatically exist in PostgreSQL 9.2? because I also need to use pg_routing


texts: Dynamically Create a Foreign Key for Relationship Classes

texts:
May be a relativity novice question but my research has turned up no direct answer.

I have a feature class of land owner parcels that I would like to preserve its schema. I would like to add fields that users could edit to determine various information about the parcel. My thought was to create a relationship class with a stand alone table that would contain the added fields. It would maintain 1:1 cardinality.

Here's the question.. Is there a way to configure this so that the foreign key is dynamically created from the origin table? Such that if a new parcel is added a new row in the table is created simultaneously?


texts: Configuring HERE maps on CartoDB Builder (config/app_config.yml)

texts:
I'm installing the CartoDB platform on productive environment. I'm working on about 2 weeks with all tutorials I can find with Google, reading the repository documentation, and all stuff that can help me to finish this job.

I'm stuck to configurate HERE basemaps, I have the app_id and app_code. I'm try to configurate by file config/app_config.yml, but without successful.
I try basically everything I could find on the web but still can't show the option on CartoDB Builder.

Where can I find some documentation?


texts: Coordinates problems when converting DXF to SHP

texts:
I had a DXF converted by a QGIS plugin called "AnotherDXF2Shape V0.9". When I overlayed the shapefile on Google Maps, the draw is not in the same position as it should have been. The distances between the draw and the details on the imagery of Google Maps are several feet.
Can anyone help me understand the problem?


texts: Upload scene layer package to ArcGIS Online and show all features

texts:
I am currently working in Arcgis Pro 1.4.1, creating scene layer packages. I successfully created a package of one layer that consists of 35,000 points. However, when opening it in a scene viewer in ArcGIS Online it will only show a few random points. Zooming in and out gives me more points and a better view but still the points are shown rather random. Some points that are in front of other points pop up later than others.

I already tried to upload a package with only 2,000 points but this gives the same effect. I also tried to upload points through a hosted webscene but than I am limited to 2,000 points.

Does anyone know how to upload all the 35,000 points to a scene viewer, showing all of the points, or at least the most prominent ones first?


texts: GDAL (Python) reduce raster by systematic sampling

texts: Is there any way to reduce a gdal raster in Python by sampling e.g. every 5th cell in both longitude and latitude directions?  The methods supported by gdal (e.g. gdalwarp, this question) don't obviously support this.  Is anyone aware of a method?


texts:
When I try to use the define projection tool to change the projection of a layer it will not work for some reason. I am not sure what to do about this. I am trying to change it do GCS_WGS_84

Removing the dash in the layer name did not help, here is the full error message:

Executing: DefineProjection CULT_ARE GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]


texts: Where is the source code of QGIS Processing algorithms?

texts: I am beginning to learn QGIS 3 up from QGIS 2.  I see many changes.  In QGIS 2 one could find the scripts for the algorithms located in the Processing Toolbox.  But here in QGIS 3 I cannot find the RandomExtract.py file.  Can someone direct me to the path?


texts:
Is it possible to enable Cross-Origin Resource Sharing (CORS) in OpenLayers?  If so, how?

I have a web map in OpenLayers, and I'm using an HTML canvas to render it as a PNG image.  My script works great, but due to the browser's same-origin policy, I can only convert my canvas to an image if all the map tiles I've drawn to it come from the same host (e.g. sub.domain.com).

To boost my map's loading speed, I pull tiles from an array of subdomains: a.domain.com, b.domain.com, etc.  This greatly improves performance, working around the browser's limit of four simultaneous connections per host, but prevents me from ever converting my canvas to an image, as the tiles come from multiple hosts.

I'm trying to use CORS to work around this problem.  I've set the appropriate Access-Control-Allow-Origin response header for the tileset I serve with PHP, but this doesn't have an effect.  I suspect this is due to not having the Origin header in the request (as shown in the example on the Wiki page).  It would seem that some OpenLayers configuration is required.  But what?  Has anyone else done this successfully?


texts: Value error when trying to load line and polygon GEOJSONs into PyCharm/Python

texts:
I'm trying to use GEOJSON files to visualise some data in Python using GeoPandas (GDAL, fiona, shapely, numpy, matplotlib are all installed). I'm having success in loading some files but not others; those that can be read appear to be pointers whilst lines and polygons cannot (this is based on the names of the files).
On using geopandas.read_file I get the below error:
Traceback (most recent call last):

  File "C:/......../PycharmProjects/Mapping/Mapping Code.py", line 41, in <module>
    gambia_map = gpd.read_file(

  File "C:\Python3.8\lib\site-packages\geopandas\io\file.py", line 127, in read_file
    gdf = GeoDataFrame.from_features(f_filt, crs=crs, columns=columns)

  File "C:\Python3.8\lib\site-packages\geopandas\geodataframe.py", line 300, in from_features
    d = {"geometry": shape(f["geometry"]) if f["geometry"] else None}

  File "C:\Python3.8\lib\site-packages\shapely\geometry\geo.py", line 44, in shape
    return MultiPolygon(ob["coordinates"], context_type='geojson')

  File "C:\Python3.8\lib\site-packages\shapely\geometry\multipolygon.py", line 64, in __init__
    self._geom, self._ndim = geos_multipolygon_from_py(polygons)

  File "C:\Python3.8\lib\site-packages\shapely\geometry\multipolygon.py", line 143, in geos_multipolygon_from_py
    geom, ndims = polygon.geos_polygon_from_py(ob[l][0], ob[l][1:])

  File "C:\Python3.8\lib\site-packages\shapely\geometry\polygon.py", line 509, in geos_polygon_from_py
    ret = geos_linearring_from_py(shell)

  File "shapely\speedups\_speedups.pyx", line 401, in shapely.speedups._speedups.geos_linearring_from_py

ValueError: GEOSGeom_createLinearRing_r returned a NULL pointer

Process finished with exit code 1```


texts: Are there any freely available WMS tile servers?

texts:
I am looking for an open public endpoint to retrieve WMS map tiles from, for use with leaflet maps

For clarity, I am not looking to set up a tile server or for tile server software.

It looks like the NASA service is now closed


texts: Error when exporting individual pages from ArcGIS DDP

texts:
I'm working with a 52 page map book using data driven pages. I have used this before and haven't had any issues but I am currently getting an error message when I try and export individual pages. I can export the first initial page just fine, but then when I move to the second page and try to export that individually, it just goes into a cycle that doesn't seem to end in it actually exporting. Once I've closed the program, the error message that pops up reads "Exception thrown in destructor. Encountered an improper argument." I've never even seen this error message before and have no idea what it means.

I did talk to tech support about it and it seems to have something to do with the aerial basemap; once that is turned off things run much more smoothly. But I do need that imagery included in my maps and so far we don't have a fix.

Has anyone else encountered this issue and/or know of a possible fix?

Note: I'm using Desktop 10.3.


texts: Understanding Clip Warning 000117 for Polygons?

texts:
This error makes little sense because 95% of the time this works without any problem, but here is the situation:

I am clipping polygon A by polygons B-H.

While most of the time this works perfectly, occasionally this happens:

All layers are from the same coordinate system, and the processing extent is the same as polygon A when I setup the clip. Polygons C-H clip just fine without any errors, so it is usually polygon B that has the issue; however, most of the time there are no issues whatsoever. The polygons do overlap each other. I have repaired the geometry, but the results are the same.

Given that it isn't a coordinate problem, a processing extent problem, or even a geometry problem, what else could it be causing this error, and what would I have to do to resolve this problem?

UPDATE:

The error happened again with my largest set of Polygon A; however, I then deleted 75% of the Polygon A set and it worked just fine. So my question now is: how limited is the Clip feature when it comes to large areas, and is there a way around this without having to divide the Polygon A shapefile into several segments?


texts: Return values from combining different fields

texts:
I'd like to classify parcels according to combinations of front setbacks and side setbacks. I have a "Fr_sb" field with values 0 (no front setback) and 1 (with front setback), and a "Si_sb" field with values 0 (no side setbacks), 1 (setback on one side of the parcel) and 2 (setbacks on both sides of the parcel). I'd like to combine those values on a third field "Setbs", as follows:

if "Fr_sb"=0 and "Si_sb"=0 return "0"
if "Fr_sb"=0 and "Si_sb"=1 return "1"
if "Fr_sb"=0 and "Si_sb"=2 return "2"
if "Fr_sb"=1 and "Si_sb"=0 return "3"
if "Fr_sb"=1 and "Si_sb"=1 return "4"
if "Fr_sb"=1 and "Si_sb"=2 return "5"

I'm still learning to code though, and have no idea how to do this using the Field Calculator.


texts: QGIS: How to add a second polygon that is exactly 5km larger in all directions?

texts:
Given is a map with an airport. Within a distance of 5km around the runway there is a so called no fly zone for drones (in Switzerland). I want to show this zone in a map.

What I would like to do is the following:

Add a polygon around the runway (that's ok, exept that I prefer to enter the exact geo-position in decimal format instead of clicking the map).
Add a second polygon that is exactly 5km from each side away
Crop the corners, because these are now more than 5km away from the runway corner. To simplify one more point at 45° would be ok; I don't need a curve.

How can this be achieved with QGIS? Or better with another tool?

I have tried using buffer as suggested in the answer and this is the result:

Why does the buffer appear distorted instead of a uniform 5km?


texts: Removing the duplicate values in QCombo Box

texts:
In my plugin, there is a combo box which lists city names for a particular state. FOr example, if i chose a particular state and city will be some 20 names in the dropdown. But the first name alone appearing in the dropdown is repeated thrice. Other names are once in the dropdown.
The fuction which i wrote for listing the city names in combo box is shown below.

 def findTheCityNames(self):
    layers = self.iface.mapCanvas().layers()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            layers = self.iface.activeLayer()
            if layer > 0:
                self.dlg.Citycombo.clear()
                ct = self.dlg.Statecombo.currentText()
                st = "name_1='"+ct+"'"
                exp = QgsExpression(st)
                fidAndVal = {}
                #QMessageBox.information(self.iface.mainWindow(), "info", st)
                layer = self.iface.activeLayer()
                iterFeat = layer.getFeatures()
                for f in iterFeat:
                    val = exp.evaluate(f)
                    if bool(val):
                        attrs = f.attributes()
                        idx = layer.fieldNameIndex('name')
                        cityName = f.attributes()[idx]
                        fidInt = f.id()
                        fidAndVal[fidInt] = cityName
                        for x in fidAndVal.keys():
                            self.dlg.Citycombo.addItem(fidAndVal[x])


texts: ArcGIS python script validation to update text based on list?

texts:
I'm having issue with using validation to update an input box. Some background:

I created a tool to update and resize text in ArcGIS. This is done by having a series of text elements and shape elements with matching names. So there is a shape element named "Title" and a text element named "Title". I can then resize and move the shape element around in my layout and then run this tool to update the text, resize and center within the shape element.

The tool has a validation script that is supposed to generate a list of these element names in box A and then based on which one the user selects, place the actual text in box B. In this way the user can see exactly which text element they have selected and then edit that text in box B before running the tool. The problem is that after you edit the text in box B, the validation resets the text to the original value.

How can I prevent this?

This is what I tried:

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""

    mxd = arcpy.mapping.MapDocument("CURRENT")

    if self.params[0].altered:
        temp=arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT",self.params[0].value)
        self.params[1].value = temp[0].text

return

So the behavior now is that the first box populates without issue but switching between values doesn't update box B with the text.


texts: QGIS Plugins - DLL load failed: %1 is not a valid Win32 application

texts:
This is my first question so please go easy.
I am trying to develop a Plugin for QGIS. I made my plugin using the "plugin builder", then compiled the .qrc and .ui files with PySide-rcc.exe and PySide.uic.py respectively. This all seems satisfactory.

However, my first problem arises when the plugin loads as the QGIS and Python (PySide) libraries are not in the same place.

ImportError: No module named PySide

I resolved this by adding a reference to the PySide path at the start of the compiled .ui module:

Import sys
Sys.path.append("C:\\Python27\\Lib\\site-packages")

Unfortunately when I run it again I get the critical error:

ImportError: DLL load failed: %1 is not a valid Win32 application.

Google searches are leading me to believe that this is caused by an incompatibility between 32bit and 64bit modules, however, I cannot easily download different versions as they are supplied under my companies IS policy.

This seems to apply to all the QGIS modules as i can not access any modules even from the python console:

>>> import sys
>>> sys.path.append("C:/Program Files (x86)/Quantum GIS Lisboa/apps/qgis/python")
>>> for path in sys.path:
print path

C:\Windows\system32
C:\Python27\Scripts
D:\HRW\VET\HrwPyLibs
D:\HRW\VET\HrwPyLibs\Xdmf
D:\HRW\VET\HrwPyLibs\Converters
C:\Python27\python27.zip
C:\Python27\DLLs
C:\Python27\lib
C:\Python27\lib\plat-win
C:\Python27\lib\lib-tk
C:\Python27
C:\Python27\lib\site-packages
C:\Python27\lib\site-packages\PIL
C:\Python27\lib\site-packages\win32
C:\Python27\lib\site-packages\win32\lib
C:\Python27\lib\site-packages\Pythonwin
C:/Program Files (x86)/Quantum GIS Lisboa/apps/qgis/python
>>> from qgis.core import QgsFeatures

Traceback (most recent call last):
File "<pyshell#5>", line 1, in <module>
from qgis.core import QgsFeatures
ImportError: DLL load failed: %1 is not a valid Win32 application.

Has anyone else experienced this problem or is there an obvious issue that i have not noticed?

Thanks,

Chris

Running: Windows 7 (64bit), Python 2.7.2 (64bit), QGIS 1.8.0 (???)


texts: Using sys.exit(0) in an ArcMap Script Tool?

texts:
I'm using a ArcMap script tool to test if a facility is on a list. If it's not on the list, I want to print a statement and exit the script. The script works just fine as a standalone script but it errors out when I turn it into a script tool.

import arcpy, sys

FacID_raw = arcpy.GetParameterAsText(0)

POT_ID_List = ["AAA", "BBB", "CCC"]

for POT in POT_Search_List:
    POT_ID_List.append(str(POT.getValue("Facility_ID")))

if Fac_ID not in POT_ID_List:
    print Fac_ID + " not in List. Notify team lead."
    sys.exit(0)

The error I get is:

Error Info:
      <type 'exceptions.SystemExit'>: 0

How do I exit out of a script tool when a condition isn't met?


texts: Is there a method to create offline geocoder using OSM xml data?

texts:
I am trying to develop an offline android map that makes use of OSM map tiles of 15 zoom level, For this purpose i used an open source library : Osmdroid
Now, i m looking into the possibility of creating an offline geocoding/ reverse geocoding that can be integrated with my application,
can i use Osm xml data for that purpose? if so , then can anyone suggest/explain how to use it to create SQlite db.. to be used wit my app
While doing this, i hope i m not violating any data usage terms or conditions?
Thanks


texts: How to efficiently pass QgsMapLayer to geopandas.read_file()?

texts:
I want to do statistics using geopandas/pysal on existing layers - all geopackage loaded in QGIS.
Now

layer.source()
>>>'/path/to/gpkg|layername=layer'
uri = sel_bystr(layerStr).source().split('|')[0]
lyr = sel_bystr(layerStr).source().split('|')[1].split(=)[1]
geopandas.read_file(uri,layer=lyr)

isn't beautiful is it?

layer.dataProvider().uri().uri()
>>>'/path/to/gpkg|layername=\'layer\' table="" sql='

is even worse. How can I do that 'properly'? Also any advice on marrying Geopandas and PyQGIS is highly appreciated.

Somewhat related:
By the looks of the desired formatting in this question the poster might has run into the same issue.


texts: Newly added layer has incorrect and seemingly unfixable data source

texts:
I am using ArcGIS 10.2.2. I added a layer to arc map and the layer loaded with the red exclamation mark next to it, indicating that it could not find the layer. I checked the data source and it was in a drive that the computer does not even have.

I have tried:

Changing the data source from the properties and right click/data/repair data source
removing the layer and re-adding it
Re-downloading and reading the layer

The layer is the 2007 25m raster of the Center for Ecology and Hydrology Land Cover Map. The 1990 raster has the same problem. My download is missing the 2000 raster. The 1km raster works fine. I downloaded it from digimap, to which my university has a subscription.


texts: region growing feature selection using arcpy

texts:
I have a polygon featureclass consisting of parcels touching boundary with each other and for which areas are computed in the attribute table (Figure a). I want to select polygons in an iterative way until a given threshold of area is reached.

In my system,  I want to select any one of the parcels (For example, the smallest one) as a first step (Figure b). As a second step I want to apply region growing approach that will select a neighboring parcel and add it to the present selection (Figure  c and d). However, when first three parcels are selected, their total area is equivalent to 1234 unit. I want my algorithm to stop when the area reaches , for example, 942 unit. For this, I want to split the latest selected parcel to equalize the threshold (figure e).

I am familiar with basic pyhton and arcpy stuff. If someone could give me the framework to implement this.

data = "C:\Users\ln88615\Documents\GIS_data\cadastral_computation\mask2.shp"
CadastreOut = "C:\\Users\\ln88615\\Documents\\GIS_data\\cadastral_computation\\"
mask_layer = "mask_layer"
arcpy.MakeFeatureLayer_management(data, mask_layer)

field = "Matr"
cursor = arcpy.SearchCursor(mask_layer)

for row in cursor:
   geometry = row.getValue("Shape")
   out = str(row.getValue(field))
   arcpy.SelectLayerByLocation_management (mask_layer, "BOUNDARY_TOUCHES",geometry)
   arcpy.CopyFeatures_management(mask_layer, CadastreOut + "Cadastre_" + out + ".shp")


texts: Is there any way to tell whether a park from GeoNames is small or big?

texts:
I have this query for parks in Alabama for example:

http://api.geonames.org/search?featureCode=PRK&username=genadinik&country=US&style=full&adminCode1=AL&maxRows=100

But it returns all sorts of parks like ballparks, playground, national forests.  Is there a way to distinguish city-block parks from stadiums and from big parks for outdoors?

Also, there is a value that appears named score like this: 34.15653610229492

What is that value?
Thanks!


texts: Geomajas: How do I add context menu on a map marker

texts:
I'm evaluation Geomajas. Apparently, there's very little info on this framework across the Internet. I've managed to create a context menu for the map (the menu that appears when you right click on the map).

I want to do the same for a map point marker (I may need to extend this to other geometry types, i.e., Polygons, Multipolygons, Lines, etc). Anyone with an idea on how to do this, please assist.

Additionally, when I add context menu, the "About Geomajas" menu item hides my first menu item. Is there a way that I can add menu items from below this menu, or make it the last item in the menu?

I'm adding context menu by adding the line:

map.getContextMenu().addChild(MyMenus.getMyContextMenu());


texts: Is it possible to export x and y coordinates from shapefile layer as text file in QGIS?

texts:
I am wondering if it is possible to export a list of coordinates of a shapefile layer in QGIS 3.18.
The attribute table only shows height information for every single point but no x and y coordinate.
When I select one point with the object indentifier x and y coordinates are displayed.


texts:
I have a simple script that I'm running as an ArcGIS tool that adds three fields and then calculates the fields based on expressions that derive information from existing fields. I'm struggling to understand why the script only works if it is run twice. The first time the script is run it adds the fields. The second time the script runs, the fields are populated. Any ideas on how I can get this code to work in a single run? Note: You can ignore the altercounty() function; it works and is more or less a rider on the script.

EDIT

Removed the try/except handlers as suggested and receive the following error:

  ExecuteError: Failed to execute. Parameters are not valid.
      ERROR 000728: Field TownshipNum does not exist within table
      Failed to execute (CalculateField).

Still working if it's run twice.

import arcpy

Section = arcpy.GetParameter(0)
County = arcpy.GetParameter(1)

def addfields():
    if arcpy.Exists(Section):
        fieldlist = [str(f.name) for f in arcpy.ListFields(Section)]
        if 'TownshipNum' not in fieldlist:
            arcpy.AddField_management(Section, "TownshipNum", "TEXT", "", "", "5", "TownshipNum", "NULLABLE",
                                      "NON_REQUIRED")
            arcpy.AddMessage("\n" + "Field: TownshipNum SUCCESSFULLY ADDED!")
    else:
        arcpy.AddMessage("\n" + "Field TownshipNum already exists.")

    if 'RangeNum' not in fieldlist:
        arcpy.AddField_management(Section, "RangeNum", "TEXT", "", "", "5", "RangeNum", "NULLABLE", "NON_REQUIRED")
        arcpy.AddMessage("\n" + "Field: RangeNum SUCCESSFULLY ADDED!")

    else:
        arcpy.AddMessage("\n" + "Field RangeNum already exists.")

    if 'SectionNum' not in fieldlist:
        arcpy.AddField_management(Section, "SectionNum", "TEXT", "", "", "4", "SectionNum", "NULLABLE",\
                                  "NON_REQUIRED")
        arcpy.AddMessage("\n" + "Field: SectionNum SUCCESSFULLY ADDED!")

    else:
        arcpy.AddMessage("\n" + "Field SectionNum already exists." + "\n")
    return

else:
    arcpy.AddMessage("Section not input.")

# Calculate PLSS numbers for each field
def calcfields():

    if arcpy.Exists(Section):
        fieldlist = [str(f.name) for f in arcpy.ListFields(Section)]
        arcpy.AddMessage(fieldlist)

        twnspexp = "!lndkey![6:8] + !lndkey![9]"
        if 'TownshipNum' in fieldlist:
            arcpy.CalculateField_management(Section, "TownshipNum", twnspexp , "PYTHON")
            arcpy.AddMessage("\n" + "TownshipNum successfully calculated!")

        else:
            arcpy.AddWarning("\n" + "WARNING: Field TownshipNum not in table.")
            pass
        rangeexp = "!lndkey![11:13] + !lndkey![-1]"
        if 'RangeNum' in fieldlist:
            arcpy.CalculateField_management(Section, "RangeNum", rangeexp, "PYTHON")
            arcpy.AddMessage("\n" + "Field RangeNum successfully calculated!")

        else:
            arcpy.AddWarning("\n" + "WARNING:Field RangeNum not in table.")
            pass
        sectexp = "!sectn![1:]"
        if 'SectionNum' in fieldlist:
            arcpy.CalculateField_management(Section, "SectionNum", sectexp , "PYTHON")
            arcpy.AddMessage("\n" + "Field SectionNum successfully calculated.")
            pass

        else:
            arcpy.AddWarning("\n" + "WARNING: Field SectionNum not in table.")

addfields()
calcfields()


texts: Problem using Select By Attribute tool in ModelBuilder

texts:
I tried to use the unique iterator value as input to select layer by attribute as expression then I got a syntax error.

Here is the attribute table:


texts: Relevant data for cost attribute to calculate route time

texts:
I have imported OSM data into PgRouting database.
I would like to use driving_distance() function to calculate route time.

How to add relevant time cost into database ?


texts: Why is my Nearest neighbour outputting higher value in a equally spaced grid?

texts:
So I got a point layer which is equally spaced 0,5x0,5. When I run the nearest neighbour and accounting for the extent so that the raster cells are centered, I get this result. I got maybe five areas with this big difference. The cell says -1,98 and the original pointdata is -1,5.

Shouldn't Nearest neighbour output more or less the same value of a equally spaced grid if the raster cell is centered? The output of the data is the same resolution as the grid, 0,5x0,5.

I am using Qgis for this with SAGAs tool.


texts: How can I get ALL addresses within a polygon, with Google Maps or another free provider?

texts:
I am hoping to find a good solution to doing this.

We would like to be able to define a polygon of latitude/longitude values and be able to have a dataset of all of the addresses within that polygon.

Using a free service for this would be ideal, like Google Maps.

The polygons can be the size of several zip codes and will change depending on the dataset attached to it.

Please let me know if you know how I can do this.


texts:
I recently had to reinstall my conda environment on Linux. I installed from my old environment.yml file, and it picked up the following gdal versions:
$ conda list|grep gdal
gdal                      3.2.1            py38hc0b2d6b_7    conda-forge
libgdal                   3.2.1                h38ff51b_7    conda-forge

These is fairly recent, but in 3.2.0, import as gdal was removed
Python bindings: old-style "import gdal" is no longer available.
Use "from osgeo import gdal" instead

I have tons of code that use import gdal, and it's just time consuming catching all instances import gdal and change them to from osgeo import gdal. Is there some way to automate this (other than using find and editing things by hand?)


texts: Finding line segment of point on linestring using PostGIS

texts:
Is there an easy way to retrieve the line segment of a point on a linestring using PostGIS?
To map a point onto a linestring, I use ST_(3D)ClosestPoint. Now I would like to find the index of the line segment on that linestring, the closestpoint is on.
Turf.js is returning that inline with its nearestpointonline function. That returns the line segment index, the distance between point and line, the distance along the linestring together with the closest point on the linestring in one call.


texts: Restricted Acces to the Attach file, ArcGIS Enterprise

texts: Can I restrict access to the attached files in my web applications in ArcGIS Enterprise (Portal) ? I want only some groups to have access to attached files.


texts: Grass, Route Finding with Elevation Limits

texts:
I've got GRASS (7.0) installed and I "play" with it on ocassion since I am far from a GIS pro. What I want to do is to find all road routes from a certain point limited by distance (path length, not point-to-point) and grade (elevation change per distance, e.g. 5% equals 5m rise per 100m traveled) profile. The purpose is to find bicycling routes suited to my training needs at a given time.

Since, again, I am not pro, base data must be freely available (OSM, USGS sources for Siskiyou County, CA, US) as must be any extra tools. I'm pretty sure that what I want to do is not extremely difficult; but I really have no idea where to begin. I'm reasonably competent at shell, SQL, Python, etc. coding if I've got a clue as to what algorithm to implement. Now, I need serious remedying my cluelessness. :)


texts:
I want to read soil type info from grib file, but dont know a variable name.
I expect that this variable contains a value between 1 and 9, corresponding to the possible different types of soil for each grid point.

I'm use NetcdfFile library but I can not find the variable.

Any suggestion?


texts: Creating transparent PNG tiles (for TMS service) from transparent TIFF

texts:
I read nice tutorial about making transparent DEM TIFF on (Creating transparent hillshade?)

My Goal is creating transparent PNG tiles from transparent GeoTIFF
for serving my own TMS map service.

I followed tutorial with SRTM1 N35E128.hgt dem file

gdaldem hillshade c:\n35e128.hgt c:\hillshade.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges
gdal_calc -A c:\hillshade.tmp.tif --outfile=c:\opacity.tif --calc="255-A"
gdalbuildvrt -separate c:\final.vrt c:\hillshade.tmp.tif c:\opacity.tif
gdal_translate -of Gtiff c:\final.vrt c:\final.tif

At last I got final.tif (please check attatched file)

----> https://mega.nz/#!2XYhiQoY!YPzrS0Z_9CGvIJIRkkSDlJ5M3biW9gvrGr7SHSNjYnU

then I execute

gdal2tiles -z 13 c:\final.tif c:\output

But output PNG files are not transparent at all.

(nodata option or -a 0,0,0,0 wasn't effective)

And then I open the tif file with glabal mapper and export TMS tiles
(with trnasparent png option checked) but the result was the same.
Just blank area of tile was transparent.

I want transparent tile like this. It's fully transparent.

How can I get transparent png tiles (TMS) ?


texts: Accessing OpenLayers on "Manage and Install Plugins"

texts: I am a new QGIS user and I am trying to add the OpenLayers plugins.  I went into 'Manage and Install Plugins', and made sure I had the 'All' tab selected.  However, whenever I try to search for the pluggin, it never shows up as on option.  In the dialog box I wrote out 'OpenLayers'.  Other plugins show up, but not that one.


texts: Google Earth Engine scale definition and export

texts:
I am trying to download Global Surface Water (Pekel et al 2016) data from Google Earth Engine. It has a specified maximum spatial resolution of 30 [presumably meters]. When I export the data specifying a scale of 30, it exports by default with a coordinate reference system of WGS 84. According to this question Landsat is not in WGS 84 and a scale shouldn't be specified, but as was pointed out this is a requirement to export data, but regardless, you need to specify a scale when exporting the image. A copy of this code is here (for reproducibility). I do this, and load up the image in R, and am told this:

 image <- raster("~/GSW_new/so_carolina_test.tif")
> image
class       : RasterLayer
dimensions  : 850, 944, 802400  (nrow, ncol, ncell)
resolution  : 0.0002694946, 0.0002694946  (x, y)
extent      : -81.10574, -80.85134, 32.14828, 32.37735  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : ~\GSW_new\so_carolina_test.tif
names       : so_carolina_test
values      : 0, 255  (min, max)

Indicating that the spatial resolution has been given in degrees (not meters...) and now I am not actually getting things out with fixed distances (only fixed degree distances). Worse, this doesn't match the degree spacing that's given by some of the imagery given from a separate source (.00025 by .00025): see the data documentation with a specific example here, but it is the stated nominal resolution outputted from:

var grabscale = grab.projection().nominalScale();

Additionally putting that downloaded file into ArcGIS yields the following (see screenshot). Note a resolution of 27 m - in northern latitudes this is distorted....

My understanding is that landsat data is at 30 m, and this dataset is projected in WGS84.

How can I download said 30 m data, and how can I download it over areas at 60 N or 70 N (Arctic) while accurately preserving the "30-m" resolution?

I need to accurately measure distances from water bodies, so I need accurate distance spacing.

Should I be doing CRS reprojection to UTM (my field site is approximately  62°38'56.60"N, 164° 5'40.35"W or can I download accurate data from GEE and do processing offline in R or ArcGIS?

Is the data, being in a WGS84 projection with equal spacing, no longer at 30-m resolution in the x-y directions?

EDIT
To clarify this question: Tyler Erickson noted that the dataset is in WGS84 projection, with nominal 30 m resolution at the equator. I am trying acquire the finets spatial resolution GSW dataset, and as noted above they provide their data at .00025 deg spacing while the "30 m" data from GEE is at a slightly coarser resolution. How can I, using GEE, obtain the .00025 degree spacing that is provided by the GSW developers? As part of the answer, it would be very helpful to include information (even links etc. for self-study) on how exactly these calculations are made so I can interpret these things more easily in the future. For example is this a basic formula concerning the distance between lines of latitude?

Pekel et al 2016: https://www.nature.com/articles/nature20584


texts: Large dataset, how to copy values between fields

texts:
I'm in this position: i have a feature class containing about 6 million entries. I made a summarize and and generated a new table. I joined the two tables and now I need to copy the values from the new table to the original table(only one field). I've created a new field(the same attributes as the field I want to copy from) and tried using field calculator to run  OldField = NewField on that. It processes for a few minutes then ArcMap stops working.

Any workarounds?

I've read that the software is 32 bit and can't process such large data-sets, but how can I achieve the wanted result.


texts: Producing multiple map PDF from Excel file with each record being two lat/long pairs?

texts:
I need to produce multiple maps from a single excel file consisting of approximately 1,500 records.  Each record contain two lat/long pairs (four points) and I need a PDF map of each record zoomed to the extents of the four points.  Can this be done using data driven pages or through the ArcPy mapping python object?  Thanks to all.

The records are part of a task related to public transport planning.  Four points, four lat/long pairs.  The four points are "pairs" (origin/destination and boarding/alighting). I'm using ArcGIS 10.1.  Please forgive any confusion.


texts: Drop decimal and decimal places without rounding in FME

texts:
Say I have 105.6546 and 5646.546356, in string or float fields, how can I drop only the decimal and values to the right of the decimal? So my final values will be 105 and 5646.

Particularly within FME.


texts:
I want to map the absolute increase/decrease of the population of a city. The amount of people that increased/decreased is represented by a size-related point. If it's negative I want to have it in blue, if it's positive I want to have it in red.

Does someone knows how I should do ? It should be like that :


texts:
I`m trying to implement postgis St_Project analog on C++ (calculation method is taken from here - http://www.movable-type.co.uk/scripts/latlong.html):

#include <cmath>

struct GeoPoint
{
    double longitude;
    double latitude;
};

static const double EARTH_RADIUS = 6378137.0;
static double degToRad(double deg) { return (deg * M_PI) / 180.0; }
static double radToDeg(double rad) { return (rad * 180.0) / M_PI; }
static GeoPoint projectGeoPoint(const GeoPoint& gpt, double distance, double bearing)
    {
        double angular = distance / EARTH_RADIUS;
        double bearingRad = degToRad(bearing);
        double latRad = degToRad(gpt.latitude);
        double lonRad = degToRad(gpt.longitude);

        double lat = std::asin(std::sin(latRad) * std::cos(angular) +
                                       std::cos(latRad) * std::sin(angular) * std::cos(bearingRad));

        double lon = lonRad + std::atan2(sin(bearingRad) * std::sin(angular) * std::cos(latRad),
                                                          std::cos(angular) - std::sin(latRad) * std::sin(lat));

        return  { radToDeg(lon), radToDeg(lat) };
    }

int main()
{
    GeoPoint gp = {0, 0};
    GeoPoint gp_new = MapsMath::projectGeoPoint(gp, 100000, 45);
    std::cout.precision(14);
    std::cout << "longitude_new - " << gp_new.longitude << " latitude_new - " << gp_new.latitude << std::endl;
}

Result is: longitude_new - 0.63523085373227 latitude_new - 0.6351918164606

I expect that postgis St_Project function will give exatcly same result:

select st_astext(st_project(st_makepoint(0,0,4326), 100000, radians(45)))

However, this query issues different result: longitude_new - 0.635231029125537 latitude_new - 0.639472334729198

What I`m doing wrong and how achieve exactly same result akin postgis St_Project?


texts: Rasterio - save multiband raster to TIFF file

texts:
I was trying to get Rasterio to save a 3 band raster to a TIFF file. For some reason, the example in the documentation only shows how to save raster with only 1 band. Now when I tried to save a 3 band raster, I get a really uninformative error message.

The code is below. I opened a raster to collect some of the metadata. Then I created a new downsampled version of the raster and was trying to save that to a new file. The raster has 3 bands.

dat = rasterio.open('image_original.tif')

with rasterio.open('image_original.tif') as dataset:
    data = dataset.read(
        out_shape=(int(dataset.height / 3.75), int(dataset.width / 3.75), dataset.count),
        resampling=Resampling.cubic
    )

with rasterio.open(
    'image_15cm.tif',
    'w',
    driver='GTiff',
    height=data.shape[0],
    width=data.shape[1],
    count=data.shape[2],
    dtype=data.dtype,
    crs=dat.crs,
    transform=dat.transform,
) as dst:
    dst.write(data, 3) //PROBLEM LINE

So this last line creates problems. I wrote dst.write(data,3) and get an error:

  ValueError: Source shape (1, 11404, 15902, 3) is inconsistent with given indexes 1

Alternatively if I use the code dst.write(data) then I get the message

  ValueError: Source shape (11404, 15902, 3) is inconsistent with given indexes 3

which still does not indicate what the problem is or how to solve it.


texts: Geoserver not displaying ArcSDE versioned table data as it is updated

texts:
Has anybody managed to get GeoServer 2.8 with the ArcSDE plugin (ArcSDE 10.0) to retrieve updates to data as soon as they are made when data is retrieved from a SDE view?

The SDE view is in Oracle 11.2.0.4.0 which is based on a versioned feature class that takes part in a replica, hence can not use the 'move edits to base' option.

If any features are added/edited/deleted in the database they are posted back to the SDE.DEFAULT version. But GeoServer does not display the correct data unless it is restarted or the catalog is reloaded.  It seems that it keeps querying the data at the point it makes its first connection.

I am aware the edits are stored in delta tables and not the base table.  GeoServer only has access it would seem to the base table.


texts: Getting latitude and longitude of lines in attribute table using PyQGIS

texts:
I have a road network in form of lines. I retreived the Lat/Long using field calculter in QGIS.
The method I used was:
long: x(centroid($geometry)), for lat y(centroid($geometry))

However, i want to getting the Lat/Long using PyQGIS. The PyQGIS script should automatically add the Lat/Long points in the attribut table of active layer. I tried below script but it is not working.
from qgis.utils import iface

feat = QgsFeature()
mc = iface.mapCanvas()
layer = mc.layer(0)
provider = layer.dataProvider()
provider.select()

while(provider.nextFeature(feat)):
    geometry = feat.geometry()
    print "X Coord %d: " %geometry.asPoint().x()
    print "Y Coord %d: " %geometry.asPoint().y()


texts: Exporting a DXF file to be geolocated in QGIS

texts:
My cave survey freeware can export its cave map in DXF format, where XYZ coordinates are orthogonal meters from the origin station (0,0,0 by default). A user asks me to export this DXF so that it geolocates itself properly in QGIS. a georeference is available for the origin station in the raw survey data. How can I pass this geodesic information into the DXF file for QGIS? I suggested my user to alter the origin station coordinates to match its WGS84 UTM coordinates, but that does not seem to work. He also tried changing these origin station coordinates to latitude/longitude degrees, to no avail.

Can a DXF be georeferenced so that it opens and geolocates itself correctly in QGIS without any further user intervention like using bending plugins or creating shapefiles?


texts: warning message when runing GRASS 7.4 in QGIS 3.0 for Mac 10.13.4

texts:
I'm trying to run r.reclass (GRASS 7.4) in QGIS 3.0 GIRONA in a Mac iOS HIGH SIERRA. Neverteheless, when I run the comand there is noresults at all but a warning message:

  WARNING: Default locale settings are missing. GRASS running with C
  locale. Starting GRASS GIS... ERROR: Path '//'
  doesn't exist Default locale not found, using UTF-8 Exiting...

  WARNING: Default locale settings are missing. GRASS running with C
  locale. Starting GRASS GIS... ERROR: Path '//'
  doesn't exist Default locale not found, using UTF-8 Exiting...

  Execution completed in 0.70 seconds Results: {'output':
  }

  Loading resulting layers The following layers were not correctly
  generated...

I followed all the recommendations such as Instal Python 3.6 and its manual re-direction into the QGIS system, as well as the setting of the frameworks recommended.

How do I fix this issue?


texts: ArcGIS Desktop from GetCapabilities to GetMap

texts:
When making a GetCapabilities request from ArcGis Desktop, how does it get the information needed for making the GetMap request and how is the GetMap request handled by their "sample server"?

For example (taken from their own resource center http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//005300000046000000):

If I use their own sample server and make the following GetCapabilities request:
http://sampleserver1.arcgisonline.com/ArcGIS/services/Specialty/ESRI_StatesCitiesRivers_USA/MapServer/WMSServer?version=1.3.0&request=GetCapabilities&service=WMS

But how does it get the information needed for making the following GetMap request?:
http:/sampleserver1.arcgisonline.com/ArcGIS/services/Specialty/ESRI_StatesCitiesRivers_USA/MapServer/WMSServer?version=1.3.0&request=GetMap&CRS=CRS:84&bbox=-178.217598,18.924782,-66.969271,71.406235&width=760&height=360&layers=0&styles=default&format=image/png

Is it only taken from the GetCapabilities response XML or am I missing something here? It's difficult for me to understand how ArcGis Desktop creates the GetMap request and receives the images when I'm only inputting the GetCapabilities request. Is there some documentation I'm missing somewhere? And is there some kind of guideline as to how their "sample server" might be handling the incoming requests?


texts: How to export a 16bit image from QGIS via command line?

texts:
I export portions of a height map from Qgis using the command line. The options are path, resolution and extents, but it exports an 8 bit image. I am getting terracing because of this. Is there a way to export a 16bit image or a plugin to overcome this?

C:\PROGRA~1\QGIS-2~1.1\bin\qgis --project Z:\Projects\Glacier\QGIS\Glacier.qgs --snapshot Z:\Projects\Glacier\Images\Glacier_32km_02.png --width 4096 --height 4096 --extent 280500,5380000,312500,5412000


texts: Disabling Start Editing command using ArcObjects?

texts:
How do I disable the Start Editing button in ArcMap when an added ArcSDE database version is the default?

IMxDocument pMxdoc = m_application.Document as IMxDocument;
IMap pMap = pMxdoc.FocusMap;
ILayer pLayer = pMap.get_Layer(0);
int count = pMap.LayerCount;
if (count == 0)
{
 MessageBox.Show("Layer not found in ArcMap");
}
 else
{
 IFeatureLayer pFlayer = pLayer as IFeatureLayer;
 IFeatureClass pFclass = pFlayer.FeatureClass;
 IFeatureDataset pFDSet = pFclass.FeatureDataset;
 IWorkspace pWorkspace = ((IDataset)pMxdoc.ActiveView.FocusMap.Layer[0]).Workspace;
 IVersion version = pWorkspace as IVersion;
 string versioname = version.VersionName;
 IEditor editor = getEditor();
 ICommandItem item = GetCommandItemByName("{59D2AFD0-9EA2-11D1-9165-0080C718DF97}");
 ICommand cmd = item as ICommand;
 string name = item.Name;
 int indexval = item.Index;
 bool value = Enabled();
 //ToolbarControl1.SetBuddyControl(MapControl1);
 //IToolbarPalette toolbarPalette = ToolbarControl1.GetItem(item.Index) as IToolbarPalette;
 if (versioname != "SDE.DEFAULT")
 {
   cmd.Enabled = value;
 }
  if (versioname == "SDE.DEFAULT")
 {
  Enabled();
  MessageBox.Show("Database in SDE Default version so that Start Editing option Disabled");
  }
 }


texts: sym.renderer.breakCount value not changing in ArcPy with ArcGIS Pro

texts:
Experiencing trouble with ArcGIS Pro 2.3.2 sym.renderer to assign symbology to a layer in a map project using only Python. I am automating map creation, and layer creation, symbolization. I am trying to assign symbology to a layer (polygons) using the GraduatedColorsRenderer. I update the renderer, classification method and assign it a field. Then I define the classBreakValues and the breakCount. I can only get one classification break even though I have defined 5 or more.

It appears that the sym.renderer.breakCount will not allow me to assign any value to it. Whenever I assign it to any integer it seems to not take it. When I print out the value of the sym.renderer.breakCount it is always set to 1 regardless of what I set it to. For the classification methods I have tried ManualInterval, EqualInterval, Quantile.

Here is my code:

for lyr in m.listLayers():
if lyr.name == "Binning Analysis":
    sym = lyr.symbology
    sym.updateRenderer('GraduatedColorsRenderer')
    sym.renderer.colorRamp = aprx.listColorRamps('Blues (Continuous)')[0]
    sym.renderer.classificationMethod = 'EqualInterval'
    classBreakValues = [5, 15, 35, 60, 1000]
    classBreakLabels = ["1 - 5", "6 - 15", "16 - 35", "36 - 60", "> 60"]
    sym.renderer.breakCount = len(classBreakValues)
    #sym.renderer.breakCount = 5
    count = 0
    for brk in sym.renderer.classBreaks:
        brk.upperBound = classBreakValues[count]
        brk.label = classBreakLabels[count]
        count += 1
    lyr.symbology = sym

I also put this code in the ArcGIS Pro Python window with a bunch of debug print statements to make sure all of my assignments were valid. Every thing seems valid. But the one thing that keeps happening in the Python window:

sym.renderer.breakCount = 5       #### Or any other number does not matter
print(sym.renderer.breakCount)
1                                 ### always outputs 1 no matter what I assign

What is wrong?


texts: What books, journals, and electronic resources are most valuable for expanding knowledge of GIS?

texts: What books, journals, and/or electronic resources have you found most valuable for expanding your knowledge in the GIS field, and why?


texts: Unable to writeRaster for signature "rasterPCA", "character"

texts:
I am currently trying to obtain two rasters that are the Principal component axes 1 and 2 for a list of bioclimatic variables that I am using. I am using the rasterPCA() function that exists in the RStoolbox in R, but I am unable to write it to a raster. Below is the code and the error that I am getting:

library(raster)
library(RStoolbox)

##Loading the variables and then creating a raster stack, followed by a PCA
bio3 <- raster("C:\\Users\\rameshv\\bio3")
bio4 <- raster("C:\\Users\\rameshv\\bio4")
bio5 <- raster("C:\\Users\\rameshv\\bio5")
bio6 <- raster("C:\\Users\\rameshv\\bio6")
bio14 <- raster("C:\\Users\\rameshv\\bio_14")

pres_stack <- stack(bio3, bio4,bio5,bio6,bio14)

pre_pca <-  rasterPCA(pres_stack, nComp = 2)

writeRaster(pre_pca, "C:\\Users\\rameshv\\PC.tif")

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘writeRaster’ for signature ‘"rasterPCA", "character"’

I even tried a modification of the above code

writeRaster(pre_pca,"C:\\Users\\rameshv\\4_PCAforR\\PC.asc", format="ascii", bylayer=T)

Getting the same error as above:
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘writeRaster’ for signature ‘"rasterPCA", "character"’

Any suggestions?


texts: Creating a Bing and Google Tile Layer for Leaflet

texts:
I want to use the Leaflet Javascript Library to display both a Google and Bing Tile Layer, though I'm unsure how to construct the url. The documentation shows an example using OSM:

var cloudmadeUrl = 'http://{s}.tile.cloudmade.com/YOUR-API-KEY/997/256/{z}/{x}/{y}.png',
cloudmade = new L.TileLayer(cloudmadeUrl, {maxZoom: 18});

But for some reason I can't duplicate this for either Bing or Google. Any ideas?

Thanks.


texts:
There is a NDVI raster file of a corn field. My main goal is to find the 5 percentile of the field's NDVI and create a histogram of the NDVI values. The percentile data can be extracted directly with a python code. Things get trickier when I need to generate a histogram.

Currently, to create the histogram following things are being done in ArcPy -

The NDVI image was converted to a point shapefile
The attribute table is exported in .dbf formate. The dbf file is pretty large as it contains around half million pixel values. The .dbf file structure is as follows -

The grid_code column is the NDVI pixel value.

The dbf python library is being used to read the table. Then, all values are appended in an empty list and histogram is generated from the list using matplotlibs library.

It takes almost 1-2 mins to append the values into the list. Is there any other way to just access the specific column and make a histogram out of the .dbf file. I have tried filename[2] and it is providing all values including the first cell "grid_code".

A part of the whole code is provided below -

Temp_List = []

#importing the dbf file and appending each NDVI values the temp list
Table = DBF.Table(Data_Table)
Table.open()
for record in Table:
    Temp_List.append(record[1])

#percentile calculation for further use in yield map
percentile5 = np.percentile(Temp_List, 5)

#properties of the histogram
bins = 100
n, bins, patches = plt.hist(Temp_List, bins, normed=1, facecolor='green', alpha=0.5)
plt.show(block = False)
time.sleep(5)
plt.close()


texts: I am using ArcGIS Desktop 10 and we have three different servers for accessing concurrent use licenses of ArcView, ArcEditor and ArcInfo. So, whenever I need to switch to other product, I have close currently using Catalog or ArcMap and open ArcGIS Administrator and then change the software product and License manager (Server IP) manually. So, is there any way to change the software product and License manager (Server IP) using script or batch file?


texts: MultiPatch Shapefile import into PostGIS Database

texts:
I have geometry of buildings in Shapefile MultiPatch format. I want to import it into PostGIS database in the same format.

Is that even possible? If yes, then how?


texts: How to count number of polygons inside another polygon layer?

texts:
I am looking for a tool in ArcGIS 10.1 QGIS or R that allows me to count the number of polygons that are inside specific areas (defined also as polygons).

My data consists of:
one shapefile with multiple polygons with different shapes and sizes, and one shapefile with the zone/areas


texts: Displaying only selected features on map in QGIS

texts:
Not sure if this option ever existed, but the new QGIS version hides a lot of things somewhere else and this seems so simple that I'm kind of convinced it was possible to do this before.
I have quite a large data set with ship tracks. I now want to analyse individual ship tracks (made from point layer with Points2One Plugin). However, it is so crowded that it is hard to look at individual tracks. I don't want to save each single track as a new layer, so I wonder if it is possible to select one track in the Attribute table and display only that one on the map/hide the others?


texts: One attribute category defines colour, another attribute category defines shape

texts:
I want to set the point colour to be defined by one attribute (Lithology), and the shape to be defined by another attribute (Qualifier). For example,

Lithology (RH) = Yellow, Qualifier (CR) = Square, resulting marker is a yellow square.

Lithology (AN) = Blue, Qualifier (BX) = Triangle, resulting marker is a blue triangle.

Lithology (RH) = Yellow, Qualifier (BX) = Triangle, resulting marker is a yellow triangle.

13 different lithologies, and 15 different qualifiers exist, and I have ~500 points.


texts: How to reinstall QGIS from scratch on Ubuntu?

texts:
System info: Ubuntu 16.04.3 LTS, kernel 4.10.0-33-generic

I was trying to update QGIS using these steps. For some reason, it broke in the process and I'm not able to start it anymore. Executing qgis gives me the following:

  /usr/bin/qgis.bin: symbol lookup error: /usr/lib/libqgis_core.so.2.14.17: undefined symbol: OGR_F_SetFieldNull

I didn't really know what to make of it, so I tried reinstalling the whole thing. Specifically running

sudo apt purge qgis* python-qgis* --auto-remove
sudo apt install qgis python-qgis qgis-plugin-grass

However, at the end of the installation process I now get the messages:

dpkg: dependency problems prevent configuration of qgis:
 qgis depends on qgis-providers (= 2.14.17+dfsg-1~xenial3); however:
  Package qgis-providers is not configured yet.

dpkg: error processing package qgis (--configure):
 dependency problems - leaving unconfigured
Setting up qgis-plugin-grass-common (2.14.17+dfsg-1~xenial3) ...
No apport report written because the error message indicates its a followup error from a previous failure.
                                                                                dpkg: dependency problems prevent configuration of qgis-provider-grass:
 qgis-provider-grass depends on qgis (= 2.14.17+dfsg-1~xenial3); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-provider-grass (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-grass:
 qgis-plugin-grass depends on qgis (= 2.14.17+dfsg-1~xenial3); however:
  Package qgis is not configured yet.
 qgis-plugin-grass depends on qgis-provider-grass (= 2.14.17+dfsg-1~xenial3); however:
  Package qgis-provider-grass is not configured yet.

No apport report written because the error message indicates its a followup error from a previous failure.
                                                                                No apport report written because MaxReports is reached already
                                                                                dpkg: error processing package qgis-plugin-grass (--configure):
 dependency problems - leaving unconfigured
Processing triggers for libc-bin (2.23-0ubuntu9) ...
Errors were encountered while processing:
 qgis-providers
 qgis
 qgis-provider-grass
 qgis-plugin-grass
E: Sub-process /usr/bin/dpkg returned an error code (1)

I then tried the same thing without the QGIS repositories, but with the same results, which is why I'm thinking there are still some broken leftovers from my previous attempts on my machine.

Is there any way to nuke my QGIS installation more thoroughly than by using sudo apt purge qgis* python-qgis* --auto-remove, so that absolutely nothing is left?


texts: Problem with SetValue method in ArcGIS 9.3

texts:
I am trying to make a python script for ArcGIS 9.3. I have a problem with the SetValue method. I try to use it in a very simple code but it still doesn't work : No error is return but my field is not populate. The code is the following :

import arcgisscripting
gp = arcgisscripting.create()

gp.workspace = r"Z:\PPR\PPR76\Tentative Macro Amelioration MNT"

try:

    input = gp.GetParameterAsText(0)
    field = "TEST"

    rows = gp.UpdateCursor(input)
    row = rows.Next()
    while row:
        expression = 3
        row.SetValue( field , expression)
        row = rows.Next()
except:
    print gp.GetMessages(2)

Thanks for your help


texts: Create mosaic like Voronoi Diagram from disjoint polygons

texts:
The illustration below shows the problem:

as in (a) I have a set of disjoint polygons, as geometries in PostGIS. I need something like (b), the  "mosaic" of this set of polygons, building it by a "influence region" criteria... Is like a Voronoi construction (illustrated by (c)): in fact, if the polygons was points, the influence regions are Voronoi.

Summarizing: I need a SQL algorithm (or some specific for PostGIS)  that generates the "mosaic" of a set of disjoint polygons.
(perhaps a loop of little ST_Buffer and ST_Difference operations)

PS: I need, like Voronoi's, that space delimitation (a squared frame in (b)) is ignored.

This problem is similar to this other about lines.

EDIT (after @FelixIP comment)

I preffer to stay in vector universe, to not lost precision (ex. using ST_DelaunayTriangles and adding and subtracting interiors by the original polygons, them adapting a dual graph solution)... Some simple and automatic packages like pprepair (assisted like QGIS topological tools are not automatic).
But raster is perhaps simpler and less CPU-consuming.

This "GRID process" illustration  is also valid as solution, assuming that it can allowing same precision and "euclidean influence region growing".

In ARCGIS exists a spatial analysis tool known as Euclidean Allocation, so, perhaps there are a PostGIS similar solution, starting with the set of polygons (classifying, rasterizing and making-back the polygons).


texts: OpenLayers on mobile device: possible to use local map tiles?

texts:
Is that possible for an OpenLayers app (running on mobile device) to load the map tiles that are deployed in the device, if the app is able to access these files (e.g., using PhoneGap)?
from this question:
Using local mbtiles in openlayers?
with code:
map.addLayer(new OpenLayers.Layer.TMS('overlayWorld',
'http://localhost:8888/2.0.0/overlay_62ee43/{z}/{x}/{y}.png', {
maxResolution: 156543.0339,
type: 'png',
layername: 'overlayWorld',
isBaseLayer: false
}));

it seems work. but it employs local web server. In the case of mobile device, does it still need a web-server?
any advice will be appreciated!


texts:
I am trying to convert my vector file (polygons) to raster by assigning a 10x10 size to the pixel cells.
However it gives me this error:

I used as "Field to use for a burn - in value" a column of the vector attributes (integer number).


texts:
So I am having a shapefile with all the world countries. And the table for this shapefile has a column named "Country Name". Each of this entry is unique, meaning, each country takes only 1 row, and each with it's own polygon represented in the shapefile and ID.

Now I have a .csv file with all the countries again. Now this .csv file has the column "Country Name" too, but this column has several rows containing the same country name, because, it also has a column called "Year", therefore each country has several row entries for each year.

How can I join my shapefile table with this .csv table? Right now, my problem is that when joining, only the first entry from the .csv file is being joined, since the shapefile table only has 1 row per each country entry.

Any suggestions how I can a union join, meaning that each country in my shapefile will be assigned every year found in the .csv file? I think it means that each polygon needs to be multiplied for every year entry from the .csv file.

Thanks

EDIT: According to ArcGIS Resources - joining and relating tables, "In all cases of 1:M joins, only the first matching record is joined and displayed in the layer's attribute table.". Which is actually the problem I have. Only the first matching record is displayed in my shapefile table.

So I believe I now need a work around? Since my .csv is having my data by year, but each year is a row (and the same country is spanning several rows for each year row), would I need to create a column for each of the years, and somehow transfer the data found in each row from my .csv file to the columns of my shapefile?


texts: Buffer polyline to kilometres in custom projection in QGIS, Python, R

texts:
I have a map containing polyline squiggles all over the world (shown below). I need to buffer these lines by 50 km either side, resulting in polygons.

I have reprojected my data to this custom projection:

+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

However, when I buffer my polylines using the Vector>Geoprocessing>Buffer option QGIS reverts to WSG84 with an error bar stating that my crs is unknown.

There are two options I am considering:
(1) i) Use ESPG 3857 instead, as suggested as option 4. here (although a commenter in this thread (top comment of the top answer) states that EPSG 3857 is not in real 'm' - so I am now unsure if this option is valid). ii) Reproject to my custom crs.

(2) Try to force the buffer to work in my projection using code in R or Python.

So my questions are:
a) Is (1) a valid method?
If so, despite the 'metre units not being real metres' comment, then I am down to use the easier method, followed by reprojecting to my custom crs.
b) Or if are there some code options (2) that allow me to use my own crs, I would love to hear those too!


texts:
I'm confused about the following basic points:

Do lat/long points have a projection or are they projection agnostic?
Do these points live in a sphere or a plane?

I understand that projections are used to transform a spheroid into a plane while attempting to minimize artifacts along desired attributes (e.g. distance).

Is a lat/long measurement agnostic to projections in the sense that it captures some sort of "absolute position" on Earth?
If it is in fact agnostic (if the point lives in the spheroid), is this the reason why we use havesine to calculate distances?
How do CRS come into play when dealing with lat/long coordinates?

I seek feedback on misconceptions and pointers to material to learn about this stuff.


texts: Using the MensurationDockWin in ArcObjects .Net

texts:
I have a DockableWindow add-in that had it's own map control.  I would like to be able to use the mensuration tools from the Image Analysis window.

I can add the standard Measure Tool with the following code, but it doesn't have desired tools. Namely Point Location and Base to Top Height.

ControlsMapMeasureTool measureTool = new ControlsMapMeasureTool();
measureTool.OnCreate(axMapControl1.Object);
axMapControl1.CurrentTool = (ITool)measureTool;

I've tried looking for how to create/use the MensurationDockWin, or even the ImageAnalysisWindowCommand, but I'm not sure how to go about utilizing Esri Dockable windows.

I've tried the following, but whenever OnCreate is call it throws an AccessViolationException [Trying to read/write protected memory]

MensurationDockWin mensurationWindow = new MensurationDockWin();
mensurationWindow.OnCreate(axMapControl1.Object);

ImageAnalysisWindowCommand window = new ImageAnalysisWindowCommand();
window.OnCreate(axMapControl1.Object);


texts: Error with publishing Sort tool on ArcGIS Server 10

texts:
I have a model that has one tool (Sort), works fine on desktop. When I try to publish get this error.

 Geoprocessing service initialization failed.
Invalid parameter: Tool - General Sort : Parameter - Field(s) : Data Type - Value Table

Turns out this data type is not supported, looking for a workaround that would allow users to choose more than one values for Sort field. Thanks, Jay


texts: How do you pair points on separate layers with similar attributes in Qgis?

texts:
I have two separate shape files with point data.  One is POI data and the other is exit location.  The poi data contains the unique exit id linking it to the exit.  How do I visually show this in qgis?  Some exits have multiple poi's some have none.

I've tried Points2one but it only links a few of them and it seems pretty limited.


texts: Why does my raster layer disappear when changing colors?

texts: When I change the colors of a raster layer (EVI), the layer disappears. Here are my color settings:


texts: Exporting Dot Density Map as Points using ArcGIS Desktop?

texts:
I have a polygon with a sum_count field. Using the Dot Density symbology in ArcMap, I was able to display this value as a series of points (seen in the image below).

Is it possible to export these points as their own feature?

The objective is to do some advanced cartography with these points.


texts: Creating buffer circle x kilometers from point using Python?

texts:
I am wanting to create a circle with an inputted radius around a point of predetermined latitude and longitude. I've been trying to use the formula found in this post which uses an algorithm documented by Ed Williams found here - but the outcome always creates a rectangle (seen below after plotted on google earth).

My code looks like:

bufferRange = input("enter radius of buffer zone")
bufferRange = float(bufferRange)
angle = 0
newCoords = []

for i in range(360):

    lonDD = abs(lonDD)
    d = bufferRange/6371
    Radian = angle * 0.017453
    lat = math.asin(math.sin(latDD)*math.cos(d)+math.cos(latDD)*math.sin(d)*math.cos(Radian))
    lon=math.fmod(lonDD - math.asin(math.sin(Radian) * math.sin(d)/math.cos(lat))+math.pi,2*math.pi)-math.pi

    if 0 < angle < 90:
        lat = latDD + lat
        lon = lonDD + lon
    elif 90 < angle < 180:
        lat = latDD - lat
        lon = lonDD + lon
    elif 180 < angle < 270:
        lat = latDD - lat
        lon = lonDD - lon
    elif 270 < angle < 360:
        lat = latDD + lat
        lon = lonDD - lon
    elif angle == 90:
        lat = latDD
        lon = lonDD + lon
    elif angle == 180:
        lat = latDD - lat
        lon = lonDD
    elif angle == 270:
        lat = latDD
        lon = lonDD - lon
    else:
        lat = latDD + lat
        lon = lonDD

    lon = -lon

    if angle < 360:
        newCoords.append((lon, lat), )
    else:
        newCoords.append((lon, lat))

    angle = angle + 1

Where latDD and lonDD are the points at the center of the circle (Ottawa for this test).

Update:

After using xunilk's code (see below) the results improved, but the buffer is still not circular. I have a feeling its how I converted my longitude to meters, but I can't seem to find a better way to convert them online.

The following test was done with a buffer radius of 100km. The distance from the centre to the 'top' and 'bottom' of the buffer is exactly 100km. The distance from the centre to the sides of the buffer is 70.5km

polygonSides = 360
lat = (latDD*111.320)*1000
lonNum = 111320*math.cos(latDD)
lon = abs(lonDD)*111320

points_list = [ (-1*((lon + np.sin(angle)*bufferRange)/111320), (lat + np.cos(angle)*bufferRange)/111320)
for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ]

Image of my new buffer:


texts: Adding updates to current QGIS version from later ones?

texts: I am currently using version 3.4 of QGIS, and I was wondering if it is possible to implement all of the awesome changes from version 3.5, 3.6, etc. into version 3.4 or do I need to keep installing the new versions if I want access to the new features?


texts: Shapefiles plots together in a GIS, but not in R

texts:
I tried plotting these two files and a set of points (dropbox folder) in QGIS and MapInfo. They align pretty good (not perfect though).

But in R, only the data from www.naturalearthdata.com can plot along with other things like the points. My 'own' shapefiles plots nicely, but always alone (not adding on top of or "allowing" points on top of) and...

My own shapefiles are rotated (clipped from a larger file where the center has north up) and I can't figure out how to save the applied rotation in MapInfo (I've tried exporting numerous ways). Maybe this is causing my issue?

library(rgdal) library(raster)
#Loading "my" shape files:
shp.WashingtonLand<-readOGR(".","WLmap2_polyline")
plot(shp.WashingtonLand)

#Loading naturalearthdata.com data for only Washington Land
shp.greenland<-readOGR(".","GRL_adm2")
ext.washington.land <- extent(-68, -58, 79.9, 81.2)
shp.clip.WashingtonLand <- crop(shp.greenland, ext.washington.land)
plot(shp.clip.WashingtonLand)

#data points
SamplePoints <- read.csv2("SamplenumberAndCoordinates",header=TRUE)
locs <- subset(SamplePoints, select = c("Latitude", "Longitude"))
coordinates(locs) <- c("Longitude", "Latitude")
plot(locs, col="red", add=T)

Hope it makes sense, it took a while getting all the things together and simplified appropriately.


texts: Simplifying adjacent polygons using PostGIS?

texts:
I encountered a problem simplifying set of polygons that are adjacent. If I simplify each polygon separately with the Douglas–Peucker algorithm (which is used by many open source tools), the resulting polygons are usually not adjacent anymore. This problem exists, for example, when simplifying borders of countries/provinces.

Does anyone have a solution for it using PostGIS?


texts: How can I create an isohyetal map in ArcGIS 10. I have the precipitation values from a number of weather stations in my area of interest. What tools should I use to create these contour lines?


texts: Georeferencing PDF maps for use in Avenza PDF maps

texts:
I have PDF maps like the one below (normally PDF, but JPEG for this) that I would like to be able to georeference to use with Avenza PDF maps.

I don't have GIS training, but have access to Arcmap 10.1 and Global Mapper (I looked at the Avenza software, but it is out of my price range).

How involved is Georeferencing the maps? I have the scale, projections and GPS points but am unsure of how to proceed.


texts:
Apparently the lowest scale visibility in ArcGIS Online is 1:1,250.  It would really be nice if a map could zoom in closer than this.

In fact, some of the maps that I published in previous years to ArcGIS Online did and still do allow me to zoom in closer than 1:1,250.

What changed and why are we limited to this relatively course scale?


texts: How do I create this map polygon outline in QGIS

texts:
I'm quite new to QGIS and want to create a map using this technique

The map of GB is made up of lines on one layer. Each of the yellow polygons must touch or overlap each other and have the outline of the coast where it meets the lines from the map layer. I would like to end up with a polygon shapefile which would form the basis of an exclusion zone for another program.

I'm looking for tutorial or guide on how this might be achieved. If this is not available then some simple tips on what tools in QGIS to research in order to self teach.


texts: Mapnik, TileCache, Apache2 - Setup&Debug

texts:
I'm trying to setup TileCache with Mapnik working on an Apache Server. Right now I don't see any tiles in my OpenLayers Framework in the webbrowser. I really ask myself where the problem is hidden.

a) Does anyone have a good debug strategie for working with these 3 components?

b) Anyone working with that and can give a little working sample setup. Maybe my confs are wrong or incomplete.

Thanks!


texts: Center map before exporting to jpeg with ArcPy

texts:
I am trying to standardize map export for some of my maps with ArcPy. In the ArcGIS interface, I would click FullExtent before exporting (File --> Export map), so that the full map is covered in the image file.

With ArcPy, I have not understood how this works. Currently, I use arcpy.mapping.ExportToJPEG to export the maps. When I then look at these images, however, they only show part of the map. How can I make sure ArcPy centers the map first such that the full extent of the maps is shown?


texts: I was wondering if Arcgis could export a real mime type KML file and if not how to do it. I have tried converting a layer into kml however it apparently creates a kmz which I then open in Google Earth and the saved it as a KML


texts: ArcGIS 10.2 Raster Calculator can't find numpy after new local Python was installed separately

texts:
I have been working in a 60 machine GIS lab recently, which had ArcGIS 10.2 installed on every machine and working properly. I was involved with the installation of some other GIS-related software that required Python 2.7 to be installed globally on each machine (in C:\Python27). We added that path to the PATH environment variable as usual. However, now running Raster Calculator in ArcMap results in an import error for numpy.

What is bizarre to me is that ArcGIS already has its own python installed inside of C:\Program Files (x86)\ArcGIS\etc..... which does indeed have numpy (I've checked). When I open ArcMap and use the python console, import numpy works fine, which I take to mean that the python console is using the inside-of-arcgis Python installation, not the new global one.

The Raster Calculator numpy import error goes away by simply pip installing numpy into the new global Python.

Can anyone confirm that when it comes to raster calculator, ArcMap looks outside of the internal ArcGIS python and tries to use the global one?

If so, is there a good reason for this?

Is there something I am missing here?


texts: Get point layer from CSV with coordinates in meters

texts:
QGIS 3.10. Project in WGS 84/UTM zone 41N. Trying to create new point layer from CSV file, where coordinates for X and Y are given in meters (MSK-66 zone 1 - one of local CRS for Russia based on Gauss-Kruger). While doing that through usual: Layer - Add new layer - Add Delimited Text layer I'm getting only the table with missed geometry. This is the example of coordinates I have.
390803.41,1535372.18
390809.75,1535399.92

Have found one method with online reprogecting original data to WGS84 (https://geobridge.ru), than opening CSV in JOSM, creating geojson, and than back to QGIS, but that's tediously. Is there a simpler solution?
UPD. After Gabriel De Luca question.
That's what I have in the dialog. Now I got geometry, but facing the problem of wrong location of the points.


texts: Using Network Analyst in ArcGIS Desktop?

texts:
I would like to open up a new settlement with fiberglass and i want to know how much route in meter I need to open up each house in this settlement.

For this purpose I will use the network analyst in ArcGIS.

The problem with this tool is that it always calculates the distance between the settlements entrance and each house.

The point is that there should be one route entering the main street in this settlement and from this main route there should be little routes to each house.

The main route should be only calcuated once and not for every house.

Does anybody have a solution?


texts: What tile/tiling strategy should I use with GeoNode/GeoServer?

texts:
I've got a wicked GeoNode instance running, loaded with some basic vector data. I want to be able to serve our local-County aerial imagery as a basemap within GeoNode. I am thoroughly overwhelmed with where to begin...so I'm going to lay out some of my details:

Countywide aerial imagery consists of 7 MrSID files, about 1+GB each, in California State Plane Coordinates.

From the basic research I've done, here's the strategy I was roughly thinking:

gdalbuildvrt to mosaic .sid's
gdal2tiles - would i only create one layer of tiles here  in order for ImageMosaic plugin (GeoServer) to consume?
gdalwarp - I'll need to warp tiles into Web Mercator sometime, right?

Like I said, I'm a little overwhelmed. Should I use the WMS of GeoServer or go with an independent WMS e.g. TileCache/TileStache or MapBox/TileMill? What kind of workaround will I need to do with my aerials and their native projection? Am I missing a step with gdal_retile and creation of a image pyramid? Thanks for reading, I hope this is a good question for the community ~Joe.


texts: Visualizing Clusters with Google Earth Engine

texts:
I am using kMeans clustering to visual 3 clusters of an NDVI image. The clustering is always correct, but lower values could be green, when I want them to be red. Is there anyway to guarantee consistent coloring?
I have tried hard-coding values for the classes and remapping them. Still getting random outputs.
I have added indices such as NDVI and NDRE as bands to an image collection of Sentinel-2 images.
NDVI and NDVI Clustered (Expected Result):

Using Python, but will happily accept JavaScript
Code:
# Initalise the interactive map, centered on some coordinates.
myMap = geemap.Map()
myMap.centerObject(TimeSeries_AOI, 12);
myMap.add_basemap('HYBRID')

imageList = sentinel2_AOI.toList(sentinel2_AOI.size())
listSize = imageList.size().getInfo()

'''
averageBand is a function that can get the average band in a given list.

@ band: [string] mandatory string that the function uses to display a specific band.

Returns ee.Image
'''
def averageBand(band):
    image = ee.Image(imageList.get(0)).select(band)

    for i in range(1, listSize):
        if (i > 10):
            break
        image = (ee.Image(imageList.get(i)).select(band)).add(image)

    return image.divide(i)

def displayCluster(band, numClusters = 3):
    image = averageBand(band)

    #myMap.addLayer(ee.Image().paint(TimeSeries_AOI, 0, 4), {}, 'region');
    myMap.addLayer(image, visualise[band], band)
    #myMap.addLayer(image, {}, "NDVI")

    #Make the training dataset.
    points = image.sample(**{
        'scale': 10,
        'numPixels': 40000,
        'geometries': True
    })

    #myMap.addLayer(points, {}, 'training', False)

    # Instantiate the clusterer and train it.
    clusterer = ee.Clusterer.wekaKMeans(numClusters).train(points)

    result = image.cluster(clusterer, "kMean")

    classValues = [0, 1, 2, 3, 4]
    classPalette = ['F7A738', 'CDDA49', '7DB249', '327C36', '1FBCD2']

    newResult = result.set('kMean_class_values', classValues)
    newResult = newResult.set('kMean_class_palette', classPalette)
    newResult.remap([0, 1, 2, 3, 4], [1, 2, 3, 4, 5], 1, 'kMean')

    #class_values = result.get('kMean_class_values').getInfo()
    #print(result.select('kMean'))

    #print(newResult.select("data_type").getInfo())

    #subset = newResult.select("kMean").eq(2).selfMask();

    #geemap.ee_export_image(newResult, filename="Cluster2.tif", scale=90, file_per_band=True)

    myMap.addLayer(newResult, {}, band + ": Clustered")

displayCluster("NDVI")
displayCluster("NDRE")
#displayCluster("SAVI")
#displayCluster("NDWI")
#displayCluster("CCCI")
#displayCluster("MCARI")
#displayCluster("SAVI2")
#displayCluster("MSAVIhyper")
#displayCluster("RedEdge1")
#displayCluster("RedEdge2")

myMap


texts: Display leaflet layerswitcher on click event

texts:
I would like to click on a sidebar/menu icon, and when I do the layerswitcher Icon should appear.

Right now the layerswitcher is always displayed.
I know how to make it disappear with the command

display: none !important;

I should establish a connection between the layerswitcher and the menu icon.

The menu works, so when I click their pops out a menu. The code of the sidebar menu in html is:

  &times;
  Building
  about
  <button class="dropdown-btn">
    bygning
  </button>

    Link 1
    Link 2
    Link 3

It is connected to a click event in js, where the code is:
function openSlideMenu(){
      document.getElementById('side-menu').style.width
      = '150px';
    }
    function closeSlideMenu(){
      document.getElementById('side-menu').style.width
      = '0';
    }

The js code of the layerswitcher is:

 var overlayMaps= "building": riskByg};
 var byg=L.control.layers(null, overlayMaps).addTo(map);

I assume it has something to do with the connection of my sidemenu in html. But right now I can only connect to a link, and I want to connect it to an icon.


texts: Do extra zeros change the coordinates for decimal degrees?

texts:
At my work we require contractors to submit work reports in decimal degrees with 6 decimal places, i.e. the location of where they paved roads or put in a road sign.

The question is, contractors have told me that when they take the coordinates of certain locations, there only is 4 decimal places. They argue that adding 2 zeros on the end will change the location.

i.e. 38.8897°, -77.0089° vs. 38.889700°, -77.008900°

I realize for this case it doesn't, but does this hold true for all cases?


texts: Saving layer identification metadata after export using QGIS

texts:
I have set some identification metadata like Identifier, Parent Identifier, Title and Type in a layer in QGIS and when I export the layer this changes have disappeared.
How can I solve this?
I need to have this information in my shapefile (in QGIS) to import to GeoServer and then use it as a GeoJson file.


texts: MapInfo 12 problems with PostGIS database

texts:
I am simply trying to connect MapInfo 12 to my existing PostGIS database but after hours of trying various options I have hit a brick wall and need help. So far I have:

1) Downloaded and installed the 32bit ODBC PostrgeSQL driver
2) Open MapInfo 12 and enabled DBMS Catalog via the Tool Manager
3) In the DBMS Catalog I created Create Catalog>ODBC from the drop down
4) Select Data Source clicked the Machine Data Source > New and then User Data Source then scrolled to the bottom of the list to find PostgreSQL ANSI/Unicode options - I picked Unicode
5) In the dialog box I filled in my PostGIS conneciton details and hit test, and it successfully connected
6) Back in MapInfo I choose EasyLoader from the Tools list and and hit the ODBC button in the dialogue, then Machine Data Source and from the list pick PostgreSQL
7) Back in the EasyLoader dialog at the bottom i click Map Catolog' which successfully creates the mapinfo schema and table in my PostGIS database.

Then this is where it goes wrong, it now seems I need to make each table mappable

1) Table>Maintenance>Make DBMS table mappable>ODBC connection and again choose my PostgreSQL from the Machine Data Source
2) The new dialog box opens and it has all of my schemas and tables avaliable
3) I select one of the tables and click Open
4) In the Make Mappable dialog box I have
Index Type: PostGIS
Index Column: the drop down list is empty and this is where I think there is a problem

If I choose Index Type: MapInfo the Index Column comes up with geom which is my geometry column but it also populates the X and Y column as though it is trying to find cooridnates to make the geometry.

My PostGIS geometries are standard Point, Lines, Polygons, MultiLineStrings and Multi Polygons

Does anyone know how to make a table mappable?

I am on MapInfo 12.02, PostgreSQL 9.3 64Bit and PostGIS 2.1

I have hundreds of tables to map, so surely there is a better way

Thanks


texts: Display symbols of .stylx file in WPF grid in ArcGISRuntime .net

texts:
I am trying to develop a ArcGIS Runtime WPF application. I have downloaded mil2525d.stylx file and I can apply these to symbols to FeatureCollectionLayer features using the below code. So, I could see the applied style for the FeatureCollectionLayer as a result on map using DictionarySymbolStyle  and DictionaryRenderer.

List<Field> polyFields = new List<Field>();
Field symbolset = new Field(FieldType.Text, "symbolset", "symbolset", 10);
Field symbolentity = new Field(FieldType.Text, "symbolentity", "symbolentity", 10);
polyFields.Add(symbolset);
polyFields.Add(symbolentity);
FeatureCollectionTable polysTable = new FeatureCollectionTable(polyFields, GeometryType.Polyline, SpatialReferences.Wgs84);
var path = System.IO.Path.Combine(Directory.GetCurrentDirectory(), "mil2525d(1)_1.2.0.stylx");
DictionarySymbolStyle symbolStyle = await DictionarySymbolStyle.OpenAsync("mil2525d",path);
DictionaryRenderer renderer = new DictionaryRenderer(symbolStyle);
polysTable.Renderer = renderer;

This worked fine.

But, now I would like to show all the style symbols(icons) present in DictionarySymbolStyle in a WPF grid something like below by iterating through DictionarySymbolStyle  object:

//<Grid>
//<style Grid.Row="0" Grid.Column="0">
//<style Grid.Row="0" Grid.Column="1">
//<style Grid.Row="1" Grid.Column="0">
//<style Grid.Row="1" Grid.Column="1">
//</Grid>

As I am new to ArcGISRuntime with WPF, I don't know how to proceed.


texts: How use OpenStreetMap in webpage via OSM URL without Leaflet, OpenLayers, etc?

texts:
I want to know how I can interact with OpenStreetMap in a webpage via (only HTML, CSS, JavaScript) without Leaflet, OpenLayers, MapBox or anything like them (in backend or frontend). I want use pure OpenStreetMap only through URL.

For example in below code it seems I didn't run Leaflet but when I went to inspect webpage by DevTools I got Leaflet is running! And zoom control panel belongs to Leaflet.

I learned and used this question "openstreetmap-embedding-map-in-webpage-like-google-maps" asked in stackoverflow already.

<iframe width="400" height="600" frameborder="0" scrolling="no" marginheight="0" marginwidth="0" src="https://www.openstreetmap.org/export/embed.html?bbox=51.2016004294,35.6079433209,51.5569619459,35.7819296463&amp;layer=mapnik" style="border: 1px solid black"></iframe>

But when I changed code to different coordinate and load webpage again I saw map load perfect but together Leaflet control! Whereas I didn't write any code in page it want run Leaflet!


texts: Finding Median, Max and Min values from Raster Layers in QGIS

texts:
I currently have 10 years worth of January raster datasets (working on just one month at the moment) of cloud fraction cover from MODIS, and am trying to create new 3 separate rasters which will give me the Mean, Maximum and Minimum values for all 10 years.

I am currently processing this in QGIS and have so far only been able to obtain the Mean layer by using the Raster Calculator by first summing all 10 layers (years) for January and subsequently dividing it by 10, to obtain a Mean value layer for 10 years.

I can't seem to find similar solutions here of how to obtain obtain the Maximum and Minimum values?


texts: Cannot get view and extent of GeoServer WMS on Openlayers 4

texts:
I'm trying to get a satellite image from GeoServer. I've specified the extent and projection that was given in the URL (not shown in the code). But when it renders, the image is nearly invisible and put to the corner, which means that the view(?) didn't consider the extent, center, and zoom that was specified. I've broken down the code be more specific. Does anybody know how to fix this problem?

var projection = new ol.projection({
  code: "EPSG:25832",
  units: "m"
});

var img_source = new ol.source.Tilewms({
  url: "http://localhost:8080/geoserver/basic/wms"
  params: { "LAYERS": "" } //this was specified in my actual code
});

var extent = [ minX, minY, maxX, maxY ]; //taken from URL and specified in code
var x_center = (extent[0] + extent[2]) / 2;
var y_center = (extent[1] + extent[3]) / 2;

var view = new ol.view({
  center: [x_center, y_center],
  zoom: 12,
});

var layer = new ol.layer.tile({
  visible: true,
  source: img_source,
  extent: extent,
  projection: projection
});

const map = new ol.map({
  layers: [ layer ],
  extent: extent,
  view: view
});


texts: Load layers for viewing in QGIS 3 from delimited text

texts:
I am trying to use part of a script to add points from a delimited text file, however, when I use it the layers are not loaded on the screen and not saved in the path indicated as .shp

import os
import qgis.core

indir = '/home/leo/Área de Trabalho/sle/'
outdir = '/home/leo/Área de Trabalho/sle/'
for root, dirs, files in os.walk(indir):
    for file in files:
        if file.endswith('.txt'):
            fullname = os.path.join(root, file).replace('\\', '/')
            filename = os.path.splitext(os.path.basename(fullname))[0]
            uri = 'file:///%s?crs=%s&delimiter=%s&xField=%s&yField=%s&decimal=%s' % (fullname, 'EPSG:4326', ',', 'field2', 'field1', '.')
            layer = QgsVectorLayer(uri, 'pontos_' + str(filename), 'delimitedtext')
            QgsVectorFileWriter.writeAsVectorFormat(layer, "pontos_" + str(file) + ".shp", "utf-8", layer.crs(), "ESRI Shapefile")
            QgsProject.instance().addMapLayer(layer)

when I do:

print(layer)
<QgsMapLayer: 'pontos_caminho_partículas_3' (delimitedtext)>

In this case, I want to read several .txt files, with coordinates and add these points as shape in QGIS. I am using QGIS 3.10.


texts: Iterate features on SpatialPolygonsDataFrame in R

texts: I have a SpatialPolygonsDataFrame object. Is there way to use apply on this object rather than using a for loop to iterate features in this object?


texts: Changing blending of contour label to multiply

texts:
I would like to label my maps like these in QGIS (image created using Photoshop):
Image of contour label being shown on the contour lines with with white halo and then the colored raster layer above it is blended with multiply

Instead, it shows like this. White halos on the contour remains even if the contours are not visible, can not be blended with the layers above

Is there any way to create the first map in QGIS?


texts:
I am trying to compare a list of lat/long to a shapefile and see which points are in or out of the polygon. I have discovered my shapefile is seemingly giving points in a UTM format. Below are some of the points in it, they should be in and around Chico, Ca / Paradise , CA area .

How do I convert these in Python to lat long?

I have used this site to help get this far, but these have negative UTMs so I am confused.

(-132880.64330000058, 206974.35539999977),
(-133023.45749999955, 207491.3945000004),
(-133151.86959999986, 207716.73259999976)
(-133028.29309999943, 206412.34610000066),
(-133026.2402999997, 206414.30939999968),
(-133023.83819999918, 206416.1897)