texts:
I am trying to use r.horizon GRASS module(version GRASS 7.8) in Python script. But I couldn't able to access the grass modules. Everytime, I am getting the following error.
File "C:\OSGEO4~1\apps\grass\grass78\etc\python\grass\script\core.py",
  line 343, in handle_errors returncode=returncode)

CalledModuleError: Module run None g.gisenv -n ended with error
Process ended with non-zero return code 3221225781. See errors in the (error) output.

The script that I used:
import os
import sys
import subprocess

grass7bin_win = r"C:\OSGeo4W64\bin\grass78.bat"
grass7bin = grass7bin_win

startcmd = grass7bin + ' --config path'

p = subprocess.Popen(startcmd, shell=True,
                 stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()

if p.returncode != 0:
 print (sys.stderr, 'ERROR: %s' % err)
 print (sys.stderr, "ERROR: Cannot find GRASS GIS 7.8 start script (%s)" % startcmd)
 sys.exit(-1)
gisbase = out.strip(b'\n\r')
gisbase = gisbase.decode("utf-8")

Set GISBASE environment variable
os.environ['GISBASE'] = gisbase
os.environ['PATH'] += os.pathsep + os.path.join(gisbase, 'extrabin')
home = os.path.expanduser("~")
os.environ['PATH'] += os.pathsep + os.path.join(home, '.grass7', 'addons', 'scripts')
gpydir = os.path.join(gisbase, "etc", "python")
sys.path.append(gpydir)
gisdb = os.path.join(os.path.expanduser("~"), "Documents\grassdata")
os.environ['GISDBASE'] = gisdb

import GRASS Python bindings (see also pygrass)
import grass.script as gscript
import grass.script.setup as gsetup

location = "nc_spm_08"
mapset   = "PERMANENT"
gsetup.init(gisbase,
        gisdb, location, mapset)

gscript.message('Current GRASS GIS 7.8 environment:')

Could anyone help me out with this?


texts:
I have a shapefile with transportation lines and another file with all bus stops along each line (with an attribute which links stops with specific lines to which they belong).

I need to calculate line-specific order numbers for all stops - 1 being first stop from which the bus departs, 2 second stop, etc.

The shapefile containing stops contains duplicates of every physical stop depending on how many bus lines are using it (each of these duplicates have different "line" attribute number).

Any suggestions for how to calculate the new "stop sequence" attribute for all stops?

Open Source tools would be preferable.


texts: Calculating mean altitude of watershed for each cell of DEM?

texts:
I have a DEM of river basin, and I need to create a raster where the cell value will be the mean elevation of watershed of this cell. I have calculated watershed area (flow accumulation x cell area). But how to switch to mean elevation?
Any ideas?


texts: QGIS road style appears blocky, I want it to appear continuous. Is there a style merge function or something similar?

texts:
If you look closely and see the main roads, they break down into blocks [][][][].

I was wondering in QGIS how can I have them appear continuous and not blocked up as shown below. This is one file and not multiple ones.


texts: Using ogr2ogr for converting tables from non-public schema to GeoJSON

texts:
I've used ogr2ogr to convert tables in PostgreSQL database to GeoJSON which works fine for the public schema. The question is how can I set the schema to non-public to access my tables? I've already using ACTIVE_SCHEMA option as follow but it gives me not existing relation error.

ogr2ogr -f "GeoJSON" out.json "PG:host=localhost dbname=test user=postgres ACTIVE_SCHEMA=line" -sql "select * from Road_a"

error:

ERROR 1: ERROR:  relation "road_a" does not exist
LINE 1: DECLARE executeSQLCursor CURSOR for select * from Road_a


texts: How to clip symmetric difference shapefiles in Google Earth Engine?

texts:
I need to clip symmetric difference two shapefiles in Google Earth Engine. How can I do this?
I tried this code, but error:
    var a2 = ee.FeatureCollection(table2)
    var a1 = ee.FeatureCollection(table)

    var split_a = a2.symmetricDifference(a1)

    Map.addLayer(split_a)


texts:
I am trying to read a .csv file where each row contains information about a polygon. The polygon geometry information is stored in a shape column such as 'longitude latitude altitude accuracy; ...'.
Here is a simplified example with a single square polygon:
structure(list(field = "parcel_a", shape = "5 50 0 0.5; 20 50 0 0.5; 5 40 0 0.5; 20 40 0 0.5"), class = "data.frame", row.names = c(NA,
-1L))

How can I read this data format with sf in R? I would also be interested to know where does this format come from.


texts: How to connect GeoJSON LineString coordinates with a shortest distance

texts:
I have the following GeoJSON LineString that I generate from some data (that I don't own). The problem is that when this line is rendered, the 2nd coordinate (with a positive longitude) of the line is connected to the 3rd coordinate (with a negative longitude) by making a line that crosses the entire globe instead of taking the shortest distance.
Is my GeoJSON wrong?

{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {
                "ident": "R336",
                "type": "high"
            },
            "geometry": {
                "type": "LineString",
                "coordinates": [
                    [
                        168.783111111111,
                        48.6756666666667
                    ],
                    [
                        175.945333333333,
                        50.5152222222222
                    ],
                    [
                        -176.674275,
                        51.871075
                    ]
                ]
            }
        }
    ]
}


texts:
I do a map using Leaflet, and I'd like to change the background color of a popup. I included the following rule to change the color of the popup and popup tip in <style>. Nevertheless, the box changes its color* but the tip (triangle) remains white.
How can I solve it?

    .mypopup.leaflet-popup-content-wrapper{
        background-color: gray;
        opacity: 0.9;
    }
    .mypopup.leaflet-popup-tip {
        color: red;
 }


texts: Land cover change maps - Different extent when adding a New data frame in layout view

texts: I would like to create a map where I show land cover change between year 1996 and 2002 in Indonesia. In the layout view, I added a new data frame to add the other year, but it seems that the extent is different in the first Layer than the extent of the"New data frame". At the very beginning,I used the Raster Project to change the spatial reference (used Eckert IV) so that they would also have a linear unit in meters, and not only angular unit, which is in degrees. However, the extent of the ''New data frame'" is in meters, but the "Layer frame"" is in degrees. What do I need to do to make sure that all my maps have the same extent in each Data frame I add?


texts:
I received a list (XML) of SRSes by request (by JavaScript):

http://gis1:8080/geoserver/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetCapabilities
then parse it and get a list of EPSG codes like "EPSG:1234"

How I can now get a name (description) of specific SRS like "Pulkovo 1942 / Gauss-Kruger zone 13" for "EPSG:28413"?

...or may be I can do it by OpenLayers API?

I asked this first on https://stackoverflow.com/q/41672346/7424390 but then I saw that specfic forum exists so just in case I duplicate here


texts: Combine overlapping rasters in gdal_calc.py, taking the maximum value from each pixel?

texts:
I want to combine 3 overlapping rasters to create 1 raster. Each have the same dimensions, resolution and CRS. I want to populate the new raster with the maximum value from each of the layers. However, when I use gdal_calc.py the output only writes two of the rasters.

Any ideas?

gdal_calc.py -A ./EXTENT_30_OS_clipped_nodata.tif -B ./EXTENT_100_OS_clipped_nodata.tif -C ./EXTENT_1000_OS_clipped_nodata.tif --calc="maximum(A,B,C)" --outfile=./EXTENT_1000_100_30_OS.tif --type=Byte --NoDataValue=-9999 --overwrite


texts:
I have made manually one layout in QGIS 3.10, I access it using PyQGIS and I want to add a new QgsLayoutItemMap. I'm unable to set scale correctly. Does I miss a step? Is it a bug? Are you able to reproduce?

layout1 = QgsProject.instance().layoutManager().layouts()[0]

qgsLayoutItemMap2 = QgsLayoutItemMap(layout1)
qgsLayoutItemMap2.setId('myNewId')
qgsLayoutItemMap2.attemptSetSceneRect(QRectF(10, 10, 180, 180))
qgsLayoutItemMap2.setFrameEnabled(True)
qgsLayoutItemMap2.setFrameStrokeColor(QColor(255, 0, 0))
layout1.addItem(qgsLayoutItemMap2)

print(qgsLayoutItemMap2.scale()) # Return 0.0
qgsLayoutItemMap2.setScale(40000.0)
print(qgsLayoutItemMap2.scale()) # Return nan

The second return e.g nan after setting scale seems anormal.
Any clue? I want to use the project and not layers to control the display in the QgsLayoutItemMap contrary to this test code (see qgsLayoutItemMap2.keepLayerSet())


texts:
I have 200 locations of species presence. Now, I want to calculate the distance of nearest point for all 200 points.

For example, if there is a cluster of points A,B,C,D,E,F,G,H, then what is the distance between A and its nearest point? Similarly for each point, what is the distance between their nearest point? For A, nearest point need not to be B. How can I perform this in QGIS 3.8?


texts: User Interface file compilation using OSGeo4W shell during qgis plugin creation

texts: I am trying to build a plugin in Qgis2.6 using plugin builder. And as a part of the same, I tried compiling the .ui and .qrc file using OSGeo4W shell. It is showing an error during compilation in OSGeo4W shell. Please help me how to compile the ui as well as qrc file using OSGeo4W shell or reply me where I am wrong from the below attached error. I am new in QGis plugin creation. The plugin contents inside the folder also attached.


texts: Converting csv files into shapefiles using python


texts: ArcGIS 10 querying Oracle 10g without SDE

texts:
I'm interested to see if anyone has successfully got this setup working, and if so how many sacrifices were required and to which deities.

I am hoping to ultimately get ArcGIS Server 10 serving map services with Oracle 10g data without involving SDE at all. Right now I'm just trying to get ArcMap 10 to show me some data. Based on ESRI's documentation I believe this should be possible but the docs don't cover the problems I'm having.

My machine has the Oracle 11gR2 client and this is reportedly compatible with 10g. All machines involved are 32-bit only, so the 64-bit specific Oracle connection problems don't concern me (yet). I can sqlplus to the DB with the same connection parameters.

My table looks like this:

create table arc_ora_test_2 (
  m_id int primary key,
  attr_1 int,
  attr_2 int,
  shape SDO_GEOMETRY
);

INSERT INTO USER_SDO_GEOM_METADATA (TABLE_NAME, COLUMN_NAME, DIMINFO, SRID)
   VALUES ('arc_ora_test_2', 'shape',
   SDO_DIM_ARRAY
     (SDO_DIM_ELEMENT('LONG', -180.0, 180.0, 0.0000005),
     SDO_DIM_ELEMENT('LAT', -90.0, 90.0, 0.0000005)),
   4326);

CREATE INDEX arc_ora_test_2_spatial_idx_1
ON arc_ora_test_2(shape)
INDEXTYPE IS MDSYS.SPATIAL_INDEX;

From what I can gather, in ArcCatalog 'Add Spatial Database Connection' is exclusively for SDE-enabled databases, so I'm not using this. For the OLE DB connection I've tested both the Microsoft OLE DB for Oracle and Oracle OLE DB providers. Both succeeed the test connection to the DB. However...

Microsoft OLE DB for Oracle connection cannot list available tables, reporting 'General Function Failure'
Oracle OLE DB connection can list tables, but when I try to render the data I am told 'An invalid SQL statement was used'

This second (Oracle provider) connection gets furthest but it seems that ArcMap is not writing SQL properly, which is odd as there are no definition queries or other filters - this is simply a 'show me everything' query which presumably has a bounding box in the where clause.

Is it possible that the existence of the 11gR2 client is tricking ArcMap into thinking it needs to write 11gR2-compatible queries? Is there any significant difference between 10g and 11gR2 for spatial queries?

Any thoughts on this dilemma would be very much appreciated.


texts: Broken Plugin: EOL while scanning string literal

texts:
I've written a python plugin, and it works perfectly on my pc (Windows 7). I've zipped the whole folder and sent it to the client. they report back with this (the name of plugin's folder is LTS):

Any ideas why this happened? As far as I know they use Windows too.


texts: Maximum distance between two related shapefiles

texts:
I have two shapefiles of points, the first one with the point's name and the second one with several points and each point with a name of the first shapefile point. I need to measure the farthest point from the second shapefile for each point in the first shape.
For example I have 7 points with the name A and 8 point with the name B. I need from those 7 point with the name A, the farthest point and the distance between A and the specific point. And from those 8 points with the name B, the farthest distance between B and farthest point.


texts:
I have a folder with ~1200 aerial photos in both jp2 and tif format, of which I need about 380. I wrote a python script to get the names of the files I need, which works fine. But when the script copies the data, 1. it goes very slowly, 2. the file size for the jps file expands from 92k kb to 200k kb or more. It also creates ovr files an the 150K kb range. Is this normal? ArcGIS 10.0 all round. Thanks.

lst_names = arcpy.SearchCursor(raster_names, "", "", name_field, name_field)

for row in lst_names:
    arcpy.CopyRaster_management(from_folder + row.name, to_folder + row.name)


texts: Raster Calculator merge all DN values to '1'

texts:
The title of the question should say it all, but to elaborate...
I have a viewshed in raster format containing a mix of DN values and I
want to polygonise it. When I convert the raster output, the processing takes ages and the vector output generates loads of features.

The goal is to manipulate the raster to merge all DN values to '1' so that when I polygonise it, processing time is vastly reduced and the resulting feature count is also reduced.

It would also be good to know how to 'delete' values < '1' before polygonising (I have made DN values of 0 transparent in the screenshot).


texts:
Say I have 2 rasters with same spatial extent and same resolution and same number of classes:

Raster 1: Ground Truth.

Raster 2: Result of some algorithm.

I want to see what % of pixels have same value. To achieve this I tried following in raster calculator

outras = Con("ras1.tif" = "ras2.tif", 1,0)

I am getting error:

  Runtime error SyntaxError: keyword can't be an expression


texts: arcpy.da.SearchCursor(fc,["SHAPE@"]) has no firstPoint object

texts:
I am having trouble getting the start and end point of a line using the code below.

import arcpy

with arcpy.da.SearchCursor("Lines",["Material","SHAPE@"]) as curs:
    for row in curs:
        startptx = row[1].firstPoint.X
        startpty = row[1].firstPoint.Y

The error I get is as followed:
AttributeError: 'NoneType' object has no attribute 'firstPoint'

I am using ArcMap 10.3.1


texts: Adding AddIn to ArcMap UI automatically after AddIn installed?

texts:
For example: I create a basic MainMenu "My Menu" AddIn by VS2010, and add "isRootMenu=true", complie and install the addin, open ArcMap, I have to open the customize dialog -> commands, find the "[ Menus ]" below the Categories box, and find "My Menu" in the Commands box and drag it onto the MainMenu.
Many users worried about that the operation in customize dialog is too difficult.
Is there a way to Add AddIn To ArcMap UI Automatically After AddIn Installed?


texts: How Can I View All Databases on a PostGIS server in QGIS

texts:
When I add a new PostGIS database connection to QGIS the dialog asks for the database name.

At times I want to add a database and can't remember the exact name or I want to add a database from a different server and I'm not sure what databases are available.

Is there a way to list all databases on a server within QGIS? I know I could look this information up in pgAdmin but I'd prefer to keep it within QGIS.


texts: Accessing file geodatabase (.gdb) rasters with GDAL?

texts:
I'm working on a tool that converts raster layers to arrays for processing with NumPy, and ideally I would like to be able to work with rasters that come packaged in a .gdb without exporting them all (especially if this requires engaging ArcGIS or ArcPy).

Is this possible with the OpenFileGDB driver?

From what I can tell this driver seems to treat raster layers the same as vector layers, which gives you access to some data about the layer but doesn't give you the ReadAsArray functionality.


texts: ERROR 001243 preventing automatic publishing of geoprocessing service

texts:
I am trying to automate the process of publishing a geoprocessing service to ArcGIS for Server 10.2.2, using a python script:

import arcpy
result = r"C:\data\ContourLineFunction.rlt"
connectionPath = r"C:\Users\sof\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\arcgis on localhost_6080 (admin).ags"
sddraft = "ContourLineFunction.sddraft"
sd = "ContourLineFunction.sd"
serviceName = "ContourLineFunction"

# create service definition draft
arcpy.CreateGPSDDraft(
    result, sddraft, serviceName, server_type="ARCGIS_SERVER",
    connection_file_path=connectionPath, copy_data_to_server=False,
    folder_name=None, summary="lorem ipsum", tags="gp",
    executionType="Synchronous", resultMapServer=False,
    showMessages="INFO", maximumRecords=1000, minInstances=1,
    maxInstances=5, maxUsageTime=100, maxWaitTime=10, maxIdleTime=180)

analyzeMessages = arcpy.mapping.AnalyzeForSD(sddraft)

The result from the AnalyzeForSD function contains multiple errors 001243:

{
  'errors':
    {
      (u'ERROR 001243: The ContourLineFunction/in_string_version parameter is missing a syntax dialog explanation in the item description', 92): [],
      (u'ERROR 001243: The ContourLineFunction/yCoordinate parameter is missing a syntax dialog explanation in the item description', 92): [],
      (u'ERROR 001242: Tool ContourLineFunction is missing item description summary', 80): [],
      (u'ERROR 001243: The ContourLineFunction/xCoordinate parameter is missing a syntax dialog explanation in the item description', 92): [],
      (u'ERROR 001243: The ContourLineFunction/radius parameter is missing a syntax dialog explanation in the item description', 92): []} ...

I know I can get rid of the errors by manually entering item descriptions in the Service Editor in ArcMap/ArcCatalog.

I have tried to find the settings in the sddraft file, but nothing looked like the right place for the values.

I have tried to find something inside the .rlt-file (geoprocessing result), but this was all binary stuff inside a ZIP file :-(

How can I automate the setting of the item descriptions?


texts: Why did curve in parcel fabric invert itself?

texts:
I'm using the Parcel Fabric from ArcGIS Desktop 10.5.1 and when I joined this parcel to the fabric, this curve inverted itself. The traverse grid is correct with a -55.00 radius. When unjoined, the parcel and the curve show up correctly. I've tried rejoining it, checking the joins of all the surrounding parcels, changing the radius to (+) to see if it would do the opposite, and regenerating parcel fabric on just this one parcel, and nothing has worked.

Why did the curve do this and how do I fix it?


texts: How to make QGIS indicate that a layer has not CRS set in layer properties info

texts: I have QGIS 3.4.5 32bits. I added to the canvas a shapefile which has no defined coordinate system, but QGIS correctly guesses is WGS84 and warns about it. However when I go to the properties of the layer (right click layer/properties/Source) the Geometry and Coordinate Reference System title indicates that the layer is in WGS84, which is not exactly true since I have not yet actively set the CRS. Is there a way QGIS gives the right information that not CRS is set for the layer? because it is not always the person who created the qgs project, the same person who is using it.


texts: Invalid Geometry when Importing Shapefile into Excel 2016 PowerMaps

texts:
I downloaded Public_School_Locations from NYC Open Data website

When I attempt to import Public_Schools_Points_2011-2012A.shp into Excel 2016 PowerMaps this is what I get

When I attempt to pick region name, I get a list

When I pick an entry, then the error message shows, i.e. invalid geometry

How do I troubleshoot? I downloaded QGIS and came across a QGIS Tutorial website, but how do I even begin to resolve this?

Please point me in the right direction.


texts: Combining rows in attribute tables in QGIS 3.4 with some strings that should be added together and others left untouched

texts:
I want to do the following with some fields in the attribute table but I have no idea how to do it:
Each row is about a parcel and its selling price, it also contains its address, cadastral ID etc.
The problem is that sometimes, one sale included several parcels, in this case the total price of the sale is included in each parcel row. As I would like to get the price/sqm of the land sale, it obviously gets some wrong result since the price is for all the parcels in one sale, not each individual parcel.
For instance, some land was sold for 1000 euros (just for the sake of explanation). The land is located 1, Olive Street. This land includes 3 parcels A, B and C with 300, 200 and 100 sqm respectively for a total of 600 sqm. In this case, the table of attributes look like this:
Parcel_A | 300 | 1000 | 1 | Olive Street

Parcel_B | 200 | 1000 | 1 | Olive Street

Parcel_C | 100 | 1000 | 1 | Olive Street

Right now if I just do 1000/300 = 3,3 euros/sqm this is of course wrong since 1000 euros is the total price of this piece of land which includes parcels A, B and C.
I would like to transform the previous table and arrive at the following result :
Parcel_A_Parcel_B_Parcel_C | 600 | 1000 | 1 | Olive Street | 1,66 (1,66 being the price/sqm of the total piece of land).
I guess I have to do something like this (just in terms of logic) : if the address (street name and street number) is the same, concatenate Parcel Name field, add area, keep just one instance (any instance since they will be similar) of area, street name and street number. And then (ideally), replace the original rows by the resulting row (since I won't need the "old" rows anymore).
Also, in case a piece of land only includes one parcel, it should of course keep it in place since there is no issue in this case.
How do I use the field calculator to perform this kind of thing (if it is even possible)?


texts: Embed a GeoTools map in an existing Java Swing application

texts:
I have a Java GUI Swing application with three frames. I want to load a map from a shapefile and draw the map on an existing application frame. For this I use GeoTools framework. I've looked at the examples and all examples have drawn the map on a separate GUI window JMapFrame.

(I've followed this link http://docs.geotools.org/latest/userguide/tutorial/quickstart/netbeans.html)

Is there a way which we can embed a map that we have created using Geotools to an existing Java Swing application?


texts: Incorporate Google Streetview into a Leaflet map

texts:
The person who asked the question at the link below way back in 2011 said, "I understand that it is possible to incorporate Google StreetView into non-Google apps".

Can anyone tell me how this is possible? I'd like to add street view functionality to my leaflet map.

Open source alternatives to Google's StreetView engine?


texts: Programmatically select a GROUP node using PyQGIS

texts: I have seen ways to programmatically select a current layer using QgsLayerTreeView().setCurrentLayer() but I cannot find a way to change the group selection. Does anyone have any tricks/hacks that I can use or point me to the relevant Qgs class that I can use?


texts: Get the geometry type of an empty ArcSDE feature class

texts:
I'm trying to use raw SQL to get the geometry type (e.g. POINT, LINESTRING, POLYGON) of an Oracle ArcSDE feature class using the ST_Geometry spatial type. If I run

SELECT SDE.ST_GeometryType(SHAPE) FROM STREETS

and the table is empty, I get zero results, even though I specified that it's a line layer at the time of creation. I've run into the same thing in PostGIS but there's a workaround to still find the geometry type. Is there an equivalent for Esri ST_Geometry?

EDIT: The feature class was created interactively in ArcCatalog and will be populated using SQL statements, i.e. INSERT INTO STREETS (SHAPE) VALUES (SDE.ST_Geometry('LINESTRING (...)', 4326)).


texts: Convert a lat, long, value table into a raster

texts:
I would like to to convert a table with latitude, long, and a value into a raster. My file looks like this:

X   Y   Z
32.85138942 -6.643054718    0
32.85138942 -7.059721474    1004
32.85138942 -7.143054826    1014
32.85138942 -7.226388177    1080
32.85138942 -7.309721529    1040
32.85138942 -7.39305488     0
32.85138942 -7.476388231    971
32.85138942 -7.559721583    0
32.85138942 -7.643054934    3092

I tried the following:

1) Using the Point to Raster tool

a. importing the CSV file into ArcGIS
b. displaying the XY data
c. saving it as a shapefile
d. using the “Point to Raster tool” to turn the shapefile into a raster

However, this only give me a series of squares that are not adjacent together, and not a continuous surface.

2) Using the XYZ2GRID2.1-Batch program; but I can't get it to work.

Is there an other way to convert this type of lat, long, values into a raster?


texts:
For testing purposes, I deployed geonetwork 3.0.1 on tomcat7 and ubuntu 14.04 with default configuration, especially h2 database.

Now I wanted to switch to postgres and followed the instructions found on several sources, e.g. here:

I edited the file /WEB-INF/config-node/srv.xml, enabling the line <import resource="../config-db/postgres.xml"/>

After that, I created a postgresql user and database 'geonetwork' in postgresql.

Then I changed the default database password in /WEB-INF/config-db/jdbc.properties with the information from the created database and database user and restarted tomcat.

I read here that after configuring this information and starting geonetwork, there would be run sql scripts automatically to initialize database structure for geonetwork.

I cannot find anything changed from my setup with the h2 database before. When I access localhost:8080/geonetwork, there are still layer information that I created in h2 database before. I would guess they would not be transferred.

How can I check, if I successfully started geonetwork with the underlying postgres database?


texts:
Using QGIS, I have two layers of points where most of the points are overlaying each other.I want to delete the points of the first layer that are overlaying with the points of the second layer and keep the ones that stand alone.

In ArcGIS, this tool is called "delete identical", where I have one input layer and one layer that is compared to that input.
Do you know something similar in QGIS? I searched a lot but couldn't find anything similar


texts: How to buffer a point in a Python Add-In by ArcMap's selection tolerance?

texts:
In this answer to another Python add-in question I suggested supplying a search_distance to the SelectLayerByLocation function in order to allow points and lines to be selected by an intersecting point, interactively clicked by the user using the onMouseDownMap method of a Python add-in Tool class.

In ArcObjects you would normally use ArcMap's selection tolerance in pixels and expand your point's envelope by that amount in map coordinates, but I can't find any equivalent way of doing this in arcpy.

Is there a way to access ArcMap's selection tolerance setting without
using ArcObjects? And if so, how can you convert this amount to map coordinates so that you can buffer the point or select within a search distance?

Or alternatively and more generally stated:

Is there a scale-dependent way to to expand/buffer a point so that
interactive selection of features with Python add-in Tools is a
similar user experience to the usual interactive selection tools in
ArcMap?


texts: Geoserver basic authentification problem

texts:
I have Geoserver 2.1.3 set up for WFS transaction. I am using basic authentification and I have set up the following:

A role(ROLE_RMP_WFS)
A user(user_rmp) which belongs to the role above
Data security where rmp..w is set to the role above and rmp..r is set to *
No service security
Catalog mode config set to mixed

I can query the layer with WFS just fine.

However, doing a insert post transaction through OpenLayers with the URL _http://user_rmp:rmppassword@kart10utv.ad.skogoglandskap.no/geoserver_ekstern/wfs/
fails with the following message:

<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport version="1.0.0"
xsi:schemaLocation="http://www.opengis.net/ows    http://kart10utv:80/geoserver_ekstern/schemas/ows/1.0.0/owsExceptionReport.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"   xmlns:ows="http://www.opengis.net/ows">
<ows:Exception exceptionCode="InvalidParameterValue">
<ows:ExceptionText>Error performing insert: Cannot access rmp_kartobjekt_flate_v with the current privileges</ows:ExceptionText>
 </ows:Exception>
</ows:ExceptionReport>

However, the object is inserted just fine running the same query in CURL:
curl -vv --header "Content-Type:text/xml" -X POST -d @feature.xml _http://user_rmp:rmppassword@kart10utv.ad.skogoglandskap.no/geoserver_ekstern/wfs

The feature.xml is the same as the POST content:
http://pastebin.com/yPwFb9g7

I do not understand why it works using CURL, but not working using a regular browser?


texts: Raster calculator con statement not working in ModelBuilder in ArcGIS 10.3

texts:
I am trying to use raster calculator using a con statement loop for 3 rasters. All rasters are 10*10. My command is

Con(("hydro_sg2.tif" == 1), 1,
Con((("hydro_sg2.tif" == 2) & ("avg_slope1.tif" < 4)), 1,
Con((("hydro_sg2.tif" == 2) & ("avg_slope1.tif" >= 4) & ("avg_slope1.tif" <= 6) & ("kfactor_2.tif" < 0.32)), 2,
Con((("hydro_sg2.tif" == 2) & ("avg_slope1.tif" >= 4) & ("avg_slope1.tif" <= 6) & ("kfactor_2.tif" >= 0.32)), 3,
Con((("hydro_sg2.tif" == 2) & ("avg_slope1.tif" > 6)), 4,
Con((("hydro_sg2.tif" == 3) & ("avg_slope1.tif" < 2)), 1,
Con((("hydro_sg2.tif" == 3) & ("avg_slope1.tif" >= 2) & ("avg_slope1.tif" <= 6) & ("kfactor_2.tif" < 0.28)), 2,
Con((("hydro_sg2.tif" == 3) & ("avg_slope1.tif" >= 2) & ("avg_slope1.tif" <= 6) & ("kfactor_2.tif" >= 0.28)), 3,
Con((("hydro_sg2.tif" == 3) & ("avg_slope1.tif" > 6)), 4,
Con((("hydro_sg2.tif" == 4) & ("avg_slope1.tif" < 2) & ("kfactor_2.tif" < 0.28)), 1,
Con((("hydro_sg2.tif" == 4) & ("avg_slope1.tif" < 2) & ("kfactor_2.tif" >= 0.28)), 2,
Con((("hydro_sg2.tif" == 4) & ("avg_slope1.tif" >= 2) & ("avg_slope1.tif" <= 4)), 3,
Con((("hydro_sg2.tif" == 4) & ("avg_slope1.tif" > 4)), 4)))))))))))))

If I run the command just in raster calculator , it works fine ... if I use it in model builder I ma getting an error

ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 17, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 255, in Con
    where_clause)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 239, in Wrapper
    ["IfThen", in_conditional_raster, in_true_raster_or_constant])
RuntimeError: ERROR 999998: Unexpected Error

.

Any suggestions ???


texts:
So I'm trying to get the area of a buffer zone with a 50 meter buffer zone in ArcMap 10.6.1. Ideally, this would be 7853.98 meters, but when I check the Shape_Area (by using Project to change the coordinates from WGS 1984 to WGS 1984 Canada Atlas), the area is 8804.219131 meters instead. That's a pretty big difference. I checked and both coordinate systems both use meters as the base unit. What else could be wrong?

Note: The buffer zone is a small part of a process where I'm creating shapefiles from select portions of the buffer. For some reason I can't edit features created using a PCS instead a GCS, so that's why I'm having to convert coordinate systems.


texts: TextSymbol not working inside Mapview (ESRI ArcGIS 4.3)

texts:
The following is my code sample:

I tried rendering the map using the Picture Marker Symbol (as used in the example) and it works pretty well. But if I use TextSymbol as mentioned in the documentation inside Graphic object, the application fails with the error saying

  TypeError: Cannot read property 'setTransform' of undefined and u is
  undefined.

    let mapProperties: any = {
        basemap: "streets",
        ground: "world-elevation"
      };

      let map: any = new Map(mapProperties);
      let constructMap = function (lat, long, popupTemplate) {
        var point = new Point({
          longitude: lat,
          latitude: long
        });
        var markerSymbol = new PictureMarkerSymbol({
          url: "https://unpkg.com/leaflet@1.0.3/dist/images/marker-icon.png",
          width: "25px",
          height: "41px"
        });

        var textSym = new TextSymbol({
          text: "Sample Text"
        });
        console.log(textSym, markerSymbol);
        var pointGraphic = new Graphic({
          geometry: point,
          symbol: markerSymbol,
          popupTemplate: popupTemplate
        });
        var a = new MapView(mapViewProperties);

        a.graphics.add(pointGraphic);


texts: How can I get a road network to run a simulation on?

texts:
I'm new to GIS and I'm asking this a bit out of desperation, but I'm trying my best not to get overwhelmed. :)

I want to download a file that contains data of a road network (of what road network is not a concern to me), visualize this road network, but then be able to run a simulation on it (so it can't be an image, i.e. a .jpeg). I'm looking for what resources/software I should use.

(Is there a way I can chat with someone who knows about GIS? I feel that would be better.)


texts: What geographic coordinate system is this?

texts:
I am looking at a text file with data for San Francisco, CA and the "xPos" and "yPos" fields are as an example:

(-2292840, 2017170), (-2292570, 2017170), ...

I've been stuck for hours. Does anyone know what coordinate system this is using?


texts:
I have never had direct access to any Esri product, but I have written code that has interacted with ArcGIS Server, with layers that have the word "MapServer" prominently in the URL and in the JSON they provide.

But what exactly is "MapServer" (not to be confused with the other MapServer)? Is it a product? Is it just a protocol, akin to WMS? Is there actually any stand alone service called MapServer, or is it (as I suspect) just a name given to the protocol used to access generated raster map images?


texts: FME How to read multiple CSV files from many folders and write to GDP tables by folder

texts:
I have 13000+ CSV files, all containing around 13000 rows. These CSV files are arranged in 146 different folders. All CSV files have the same schema with ; as separator. These CSV files contain about 175 million rows of data combined. The folders are named as seen in the picture below (for example ...\CSV_folders\5785xxx)

CSV files are named like these:

I already tried to copy all CSV files to same directory and run them to one GDB with FME. This went well... but resulted too heavy ESRI GDB that takes forever to search through. (175 million rows, around 20GB)

Now I would need to make things a bit differently. I would like to work with FME to create one ESRI Geodatabase. It should contain tables named like these folders containing CSV files. A table in the GDB should have all the CSV files loaded in the table, that are inside the original folder. This would split my 175 million rows to 146 tables, making it a bit faster to search through with Arcmap.

So my problem is, how can I make FME read CSV files from all of the folders and write a GDB with tables named like the folder, each table containing the same information the CSV files have in the original folder combined? I'm not familiar with PostGIS etc. (and it would not work with ArcGIS), so this needs to be done with some kind of ESRI GDB workaround..


texts: Using ArcPy Result object associated with Currently Running Script

texts:
We are at the end of a Python script using ArcPy.  We can crawl around ArcPy to get information, but is there any way to get a result object for our currently running script?

My motivation is to use existing code to process the same result object we might get from executing a separate script tool earlier in the script.

To clarify, I want this object:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Result/000v000000n7000000/

I was wishing I could get a Result object for the currently running process.  Here are my notes for accessing ArcPy instead:

#Result.inputCount and Result.outputCount
outParameters = []
inParameters = []
for parameter in arcpy.GetParameterInfo:
    if parameter.direction == "Input":  #unsure of exact string
        inParameters.append(parameter)
    else:
        outParameters.append(parameter

inputCount = len(inParameters)
outputCount = len(outParameters)

#Result.MaxSeverity
maxSeverity = arcpy.GetMaxSeverity()

#Result.MessageCount
messageCount = arcpy.GetMessageCount()


texts: Snapping function doesn't work with complex polygons in QGIS?

texts:
I need to snap out parts from a polygon. I am using the "split part" to mark the area I need to erase:

In order of snapping an area, I drew a line inside the polygon. Sometimes this method does not work and I receive an error message like this:

I tried to change the snap settings, but it did not help:

Can anyone help me?


texts: Reclassify raster conditional on attribute values in R

texts:
My raster has nonsensical values (based on combination of multiple rasters in previous steps), with attributes maintained (again, due to combination of rasters in previous steps), see below for sample attribute table

I want to reclassify the Value column conditional on the 01_SHP_FINALNLCD attribute i.e. specifically, I want to reclassify all pixels with 01_SHP_FINALNLCD= 52 to 1, 01_SHP_FINALNLCD =41 to 2 etc.

I tried the following in R

How can I do this in R or arcpy? Here is the

trial<-raster[raster$01_SHP_FINALNLCD== 52 | raster$01_SHP_FINALNLCD==41]
Error in .local(x, ...) : invalid layer names

Here is the .rds file uploaded on request by @Spacedman.


texts: shapefile terminology for key-value pair

texts:
In an .osm file, there are <TAG ... /> elements like

<way visible="true" id="-4409" timestamp="1970-01-01T00:00:00Z" version="1">
        <nd ref="-21845"/>
        <tag k="fcc" v="A40"/>
        <tag k="len_miles" v="0.037850000000000002"/>
        <tag k="et_id" v="4407"/>

The "k" and "v" in these refer to key - value pairs. What are the equivalent things called in shapefiles?
Are they called (attribute, value) pairs? I see they're visible in the Attribute Table in QGIS and the k's correspond to columns and the v's are what is in the cells of that table.


texts: Add MicroStation dgn to ArcGIS/ArcCatalog 10.1

texts:
I'm an iOS programmer but new to using ArcGIS 10.1. I'm having a problem loading my MicroStation dgn file to ArcCatalog. I've connected the folder containing my MicroStation V8i roadway design dgn to the Catalog Tree and tried to define a spatial reference as described in the ArcGIS 10.1 Resource Center. When previewed, the drawing looks like a small blob that I can't zoom or view. I'm stuck to determine what the problem is. I'd think that the drawing should display in full as read only. If anyone as any suggestions I'd greatly appreciate some direction here.

Overall, I hope to overlay the dgn (drawing) over one of ArcGIS' basemaps.

I can supply any information about my dgn or setup in ArcCatalog upon request. Thanks in advance for your time.

*Here's a screenshot of how the preview of the drawing looks.

(https://www.dropbox.com/lightbox/home/Photos/Sample%20Album)

I think the blob on the right hand side is the drawing but I can't be sure what's all there because I can't zoom in to see it better.


texts: ArcMap Imagery Service Not Displaying Correctly

texts:
I am using a NAIP service to display aerial imagery. It appears as if only one tile is ever being pulled. I have played with mosaic properties of the service with no success. The tile changes every time I pan the screen, and typically covers one quarter of the screen.

Here is the service I am trying to use:

  Server: https://naip.arcgis.com/arcgis/rest/services

My co-worker does not encounter this working on Arc GIS Pro (go figure).
I am using ArcMap 10.5.1.


texts: Drawing lines between points in a shapefile that are related to each other

texts: This is a more general question, wondering if it is possible. I want to make a map of all the organizations in my city and which ones partner together. I imagine I would make a shapefile that is a list of all the organizations, and then in in the attribute table, stretch that same list across the fields, making each field the name of one of the organizations. And then for each organization that partners with another one, I could indicate that with a value. (Just to clarify, lets say organization A partners with organizations B and D. Looking along the row for the attributes of A, there would be a value in field B and D). But is it possible to display those connections on a map with lines between the organizations that partner together? So that on the map there would be a line between A and B, and another between A and D? Or is there a different way I should go about this completely?


texts: How to delete a row from attribute table on QGIS1.8.0

texts:
Is it possible to a delete row from attribute table!?
I can delete from selecting feature.

I encountered invisible polygon. I could not select feature.
Thus, I would like to know how to delete from row.

If possible,Let me know!!


texts: How to combine while and if for iteration?

texts:
I want to calculated sum of the 'Inverse_DIST' of same 'IN_FID'

I use the arcpy.da.updatecursor of python code as follows

import arcpy
arcpy.env.workspace = "D:\Auto_Logit\Neartable.gdb"
nc = max(row[0] for row in arcpy.da.SearchCursor("D:\Auto_Logit\Neartable.gdb\SP_15km_sewage", ['IN_FID']))

i = 0
j = 1
sum = 0

with arcpy.da.UpdateCursor("SP_15km_sewage",("OBJECTID","IN_FID","Inverse_DIST","Sum_Inverse_DIST")) as cursor :
    for row in cursor :
        while row[0] == j and row[1] == i :
            j = j + 1
            sum += row[2]
            if row[1] == i :
                row[3] = sum
            elif row[1] != i :
                i = i + 1
            cursor.updateRow(row)

As a result, as follows my python code not working iteration

How to doing iteration working from my code??


texts: Setting threshold for unions in PostGIS?

texts:
I have a shapefile of buffers. What I'm trying to do is to merge the buffers that intersect, but not just if any part of the buffers intersect. I would like to set a threshold, for let's say if 50% of more of two buffers intersect then join them together.

[

For example, in the representation above, I would clump the buffer zones in the bottom left, but the two clusters of two buffer zones would be kept independent.

Is there any way in PostGIS I can do this?


texts: Polygons of countries including territorial waters

texts:
Where can I get polygons (multipolygons) of all countries of the world including the territorial waters of the countries?

OpenStreetmap displays these areas but I can't find an extract of this subset of the OSM data.

is there a subset of OSM data available, just containing the territories (incl. water territories) of a country as multipolygons ?
or, can I find this set of polygons somewhere else?


texts: Installing pygame or Tkinter for QGIS Python?

texts:
I have downloaded pygame and installed it. Now i want to write scripts and run them through the QGIS console that implement methods of pygame but when the console reads "import pygame" i get a

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:/Users/stelios/Desktop/project/x.py", line 11, in <module>
    import pygame
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\utils.py", line 460, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named pygame

Any ideas?

I suppose i have to set something like a path through QGIS options but i dont know exactly.


texts: How to combine multiple line features into a single feature?

texts: I'm trying to combine line features (trail sections) into a single features.  How can this be done in ArcMap?


texts: Torque looks like a brilliant tool. However, can I use the same tool for chloropleth maps for visualizing monthly/annual changes of GDP, health expenditure, population growth etc. over the past 10 years? Will it work with bubble charts (changing the size of bubbles rather than just mapping singular events?).


texts:
I have a Leaflet map where on click I send GetFeatureInfo request to GeoServer, but the content of the HTML it returns is empty. I suppose I am computing parameters that I pass incorrectly, this is my request and how I compute parameters:

http://localhost/geoserver/cite/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=true&QUERY_LAYERS=cite:tbilisi&STYLES=&LAYERS=cite:tbilisi&exceptions=application/vnd.ogc.se_inimage&INFO_FORMAT=text/html&FEATURE_COUNT=50&X=279&Y=212&SRS=EPSG:32638&WIDTH=886&HEIGHT=811&BBOX=44.77778434753419,41.703966325182236,44.815807342529304,41.72994572390489

Code
 L.tileLayer
      .wms("/geoserver/cite/wms", {
        layers: "cite:tbilisi",
        format: "image/png",
        transparent: true,
        srs: "EPSG:32638",

        maxZoom: 20,
        version: "1.1.0",
        zIndex: 444,
      })
      .addTo(that.mymap);

    that.mymap.on("click", (e) => {
      function Identify(e) {
        var BBOX = that.mymap.getBounds().toBBoxString();
        var WIDTH = that.mymap.getSize().x;
        var HEIGHT = that.mymap.getSize().y;
        var X = that.mymap.layerPointToContainerPoint(e.layerPoint).x;
        var Y = that.mymap.layerPointToContainerPoint(e.layerPoint).y;
        let URL = `/geoserver/cite/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=true&QUERY_LAYERS=cite:tbilisi&STYLES=&LAYERS=cite:tbilisi&exceptions=application/vnd.ogc.se_inimage&INFO_FORMAT=text/html&FEATURE_COUNT=50&X=${X}&Y=${Y}&SRS=EPSG:32638&WIDTH=${WIDTH}&HEIGHT=${HEIGHT}&BBOX=${BBOX}`;
        $.ajax({
          url: URL,
          datatype: "html",
          type: "GET",
          success: function (data) {
            L.popup().setLatLng(e.latlng).setContent(data).openOn(that.mymap);
          },
        });
      }
      Identify(e);
    });

Can someone spot what I am doing wrong?
BTW a request which does generate output on the same area clicked through GeoServer Layer Preview looks like this:

http://10.28.7.8:8080/geoserver/cite/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=true&QUERY_LAYERS=cite:tbilisi&STYLES=&LAYERS=cite:tbilisi&exceptions=application/vnd.ogc.se_inimage&INFO_FORMAT=text/html&FEATURE_COUNT=50&X=50&Y=50&SRS=EPSG:32638&WIDTH=101&HEIGHT=101&BBOX=482441.61821363633,4618994.366020213,482562.1103628997,4619114.8581694765

You see at least BBOX is different.


texts: Securing one single WMS layer on GeoServer

texts:
My aim is to secure one specific layer (ex. topp:states) with username and password.

I tried a couple of combinations in the layer.properties file.

*.*.r=*
*.*.w=GROUP_ADMIN,ADMIN
topp.states.r=ROLE_WS
mode=HIDE

This somehow secures the layer topp:states but I can't access it using the role ROLE_WS (which has the user I want to give access to the layer to).


texts: What is the syntax for a Maximum operator in ArcGIS

texts:
The requested maximum operator is a function that produces a new raster with the highest values on a cell-by-cell basis from an array of coincident input rasters. It is the kind of Maximum Support classification performed by the 'HARDEN' tool of the IDRISI software.

The raster calculator spatial analyst tool should be used.

I tried:

max("raster_A","Raster_B","Raster_C","Raster_D")

But got the following error:

'exceptions.ValueError': The truth value of the raster is ambiguous


texts: Selecting multiple rows in feature based on contents of list using PyQGIS

texts:
I am trying to select rows in a feature layer in QGIS based on the values in a list. I am using PyQGIS in the python console
If I just specify one item to find in a column called "gps_code", the feature(s) are selected
I would like to pass a list of items to achieve a selection of features based on that list and have tried to use a list comprehension for this but I am not getting the desired result

layer.selectByExpression('"gps_code" = \'DAAG\'')
iface.showAttributeTable(layer)

list_of_sites = ['DAAG', 'ENTC', 'ENBR']
layer.selectByExpression('"gps_code" = \'s\' for s in list_of_sites')
iface.showAttributeTable(layer)

What am I doing wrong?


texts: How to transform Google/OpenStreetMap Co-ordinates to ESPG:3857 Co-ordinates?

texts:
I have an old map successfully layered on OpenStreetMap on QGIS and have been able to add delineated text points from a CSV file using WGS 84 / Pseudo-Mercator ESPG:3857 co-ordinates.
I have tried http://epsg.io/ set on the co-ordinate system but it isn't very easy to search for particular small road locations.
As such, I would like to use an online search tool such as Google maps or OpenStreetMaps to identify particular locations. Is there any easy way to convert these to the ESPG:3857 co-ordinate system for use in CVS files and QGIS?
For example, London's Victoria Station in ESPG:3857 is roughly -16321.992778, 6709299.869156 . On google maps, it is roughly: 51.4952074,-0.1463113 .
(From a very new user)


texts:
i have 2 polygons one is building footprints and another is boundary polygons. i want to calculate the count of building footprint which falling inside the boundary polygon. in order to avoid duplication of the building footprint count, i used 'HAVE_THEIR_CENTER_IN' option in spatial join. but the problem is that i couldn't count the footprints which fall inside the boundary polygon.

how can i count the building footprints which falling inside the boundary polygon?


texts: Packaging third-party libraries with Python Toolbox of ArcPy?

texts:
I wrote a Python Toolbox that uses libraries not included in ESRI's python (PIL, piexif).

But I can't run the toolbox on another computer that doesn't have those libraries. Since I don't want to install those libraries for every client that will be using the toolbox, is there a way to package the libraries with the toolbox to avoid any missing dependencies?


texts: Add Field disabled on Excel file in ArcMap?

texts:
I am a new user to ArcGIS Desktop 10.7. I imported data table from Excel 97-2003 format in the Table of Content. Then with XY recognition the layer was created with the defined coordinate system. The layer shows all the attributes but in the attribute table "Add Field" option does not show up. Also 'delete field', 'Calculate Geometry' and 2 other options don't activate. I tried to export the layer file to a shapefile and after the shapefile is created it does not show any attribute table.

Could anyone suggest a way to address this issue?

All I need to do is add field.


texts: Over calculation of area- How do I calculate correct acreage when there are stacked parcels?

texts:
I am currently interning for a county in the United States. I am working on the data and mapping needed to fulfill the Open Space Activity of CRS. Stacked parcels are causing me a massive headache.

As an example, take a condo community we have here. It is represented by a parcel for the land/condo association itself, and then individual parcels for each unit are stacked on top this. This is causing the acreage for such areas to be doubled, tripled, quadrupled, etc. Any ideas on how to get around this problem?

Hopefully, these pictures provide a bit more clarity of my problem. In the drawn image, the blue numbers would be parcel IDs. And represent yet another issue of smaller stacked parcels within a large parcel representing the manufactured home park or condo community boundaries.

I have also attached an image of a problem area from the actual data set. It is only a sample of a much larger data set.

I know from dissolving the feature without basing on a field (producing a single polygon) that the area should equal 26.3 acres. The stacked parcels cause the area calculation to equal about 210 acres.

FYI, if it helps, for the Open Space CRS credit I need a table with a unique parcel ID, parcel owner, parcel acreage, parcel acreage within the special flood hazard area, and open space acreage of eligible parcel.

Open space credit is determined by dividing the acreage of eligible open space by the acreage of the regulatory floodplain and multiplying by 1450. I hope this all makes sense.

https://imgur.com/a/KpJHww6


texts: GeoPandas: plot two layers but only to the extent of the smaller one

texts:
I am using GeoPandas and I want to plot two layers. I want the map to have the extent of the smaller layer. With this code below, the map's size automatically fits the bigger layer.

import geopandas as gp
import matplotlib.pyplot as plt
from shapely.geometry import LineString
from shapely.geometry import MultiPoint

thesmallpoints=gp.GeoDataFrame([[MultiPoint([(0, 0), (1, 1), (1,2), (2,2)])]],columns=['geometry'])
thelargeline=gp.GeoDataFrame([[LineString([(0, 0), (1, 1), (1,2), (20,20)])]],columns=['geometry'])

base=thesmallpoints.plot(marker='o', color='blue', markersize=20)
thelargeline.plot(ax=base, color='green')
plt.show()

How can this be achieved?


texts:
When i open file of QGIS 2.6.1 with 32 layers whose tables all come from the same database, i see that by just opening the file, QGIS already sends 4100+ "SELECT..."-questions to the database, and around 2000 "DESCRIBE..."-questions.

The interaction with QGIS becomes very slow this way.
Any idea which settings might cause this or is this by defaut?

I have pasted the SERVERLOGS on pastebin and what i see is a lot of queries like:

SELECT srid FROM geometry_columns WHERE f_table_name = 'vegetation'

http://pastebin.com/pKSgK99y

I do not have any of these columns, no "srid", no "geometry_columns", no "f_table_name" defined in my database. Should these be in database? (the qgis-database works, but very slow)


texts: How to create a polygon shapefile of most recent polygons in ArcGIS Pro?

texts: I have a shapefile that contains many overlapping polygons, each with an attribute for the year associated with the polygon (shown below). I would like to create an unstacked (no overlaps) polygon layer with polygons representing the most recent year.


texts:
I have zero experience with python but am trying to search for broken hyperlinks in a large gdb (temporarily converted to a shapefile)-

Used code from a previous question (Using ArcPy to verify paths of Hyperlinks?) and adapted it to fit my problem (didn't know if I could ask a new question on that post...)

import arcpy
import os

fc = r"G:\GISDATA\misc\project_arcpy.shp"
fields = ["PROJECT_DW"]

with arcpy.da.SearchCursor (fc,fields) as cursor:
for row in cursor:
    Roll = str(row[0])
    if os.path.exists(Roll):
        pass
    else:
        f = open(r"C:\Users\GIS\Documents\brokenlinks.txt","a")
        f.write(Roll + os.linesep)
        f.close()

Instead of writing broken hyperlinks to my text file, it writes all they hyperlinks to the file.


texts: Attribute data disappearing after converting points to line

texts:
I am currently working on a project strictly in Modelbuilder. I am trying to teach myself how to use it. I don't yet know how to use Python and am still in the early learning stage.

My question has to do with converting points to lines. I have an XY .csv file with latitude and longitude points. I am trying to convert those points to a line, but when I do so using the 'points to lines' tool, the line ends up with no attribute data in the table. All I get is one row with 'polyline M' under shape and no other values.

Here is my current workflow:

The reason that I am trying to keep the attributes is that once I create a line, I am going to measure the distance between each subsequent point on the line.


texts: MODIS data download issues with runGdal in R

texts:
I'm trying to download some MODIS data using the following function, which errors like so:
ndvi.evi <- runGdal("MOD13Q1", collection = "006", tileH = 20, tileV = 11,
                   begin = start, end = end, SDSstring = "111000000010",
                   job = "NDVI.EVI", outProj="asIn")

Error in runGdal("MOD13Q1", collection = "006", tileH = 20, tileV = 11,  :
  GDAL not installed or configured, read in '?MODISoptions' for help

When I run MODISoptions() I get the following:
> MODISoptions()

STORAGE:
_______________
localArcPath : /private/var/folders/90/qjm_6mxd1wg1f971j9t5n86c0000gp/T/Rtmpa6n6oR/MODIS_ARC
outDirPath   : /Users/u1418987/Documents

DOWNLOAD:
_______________
MODISserverOrder : LAADS, LPDAAC
dlmethod         : auto
stubbornness     : high
wait             : 0.5
quiet            : FALSE

PROCESSING:
_______________
GDAL           : Not available. Use 'MODIS:::checkTools('GDAL')' for more information!
MRT            : Enabled
pixelSize      : asIn
outProj        : asIn
resamplingType : NN
dataFormat     : GTiff
cellchunk      : 1

However when I run MODIS:::checkTools('GDAL') it seems that it can find GDAL:
> MODIS:::checkTools('GDAL')
Checking availability of GDAL:
  OK, GDAL 2.4.2 found!

I have looked at the following forum post and couldn't get this to work either, even though I successfully followed it:
GDAL not available in MODIS R package?
> Sys.getenv("PATH")
[1] "/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/Applications/RStudio.app/Contents/MacOS/postback:/Library/Frameworks/GDAL.framework/Versions/2.4/Programs"

Any advice?
P.S. When I run MODIS:::checkTools() I get the following:
> MODIS:::checkTools()
Checking availability of MRT:
  'MRT_HOME' not set/found! MRT is NOT enabled! See: 'https://lpdaac.usgs.gov/tools/modis_reprojection_tool'
Checking availability of GDAL:
  OK, GDAL 2.4.2 found!
sh: wget: command not found


texts: Clip polygons with other polygons based on matching attributes

texts:
I have two shapefiles, one which contains event locations and the the other of which represents regions where events take place.

Events polygons with year and region_name attributes
Regions polygons, which also have year and region_name attributes

The Event polygons sometimes extend beyond their proper Region.  Regions sometimes change boundaries, which is why they also have a year attribute.  I would like to clip each Event polygon using the appropriate Region polygon, i.e. the one with matching year and region_name attributes.  Here is an image I made to illustrate the basic idea:

Some solutions I've tried but which have been unsatisfactory:

Doing it manually (selecting each event and region by hand and clipping them), which seems unnecessarily tedious
Using Select By Location to isolate the Events in each Region, which still requires working with each selection manually
Using Intersect, which 1) leaves lots of event "remnants" outside their proper regions and 2) "oversplits" by dividing all Event polygons by all Region polygons, regardless of year

Note: I use both QGIS 1.8 and ArcMap 10.1.  I can't install extensions for ArcGIS (I don't have administrative access) but I can use Python scripts.


texts: Finding points near other points using ArcGIS Server 10.2 rest API?

texts:
I'm using ArcGIS server 10.2. I published my point layer in my admin. I'm going to get json file from that points which are within a fixed distance with any  points that I want.

Any idea?


texts: Modifying number of features in Merge tool each simulation run of ModelBuilder?

texts:
I am relatively new to ModelBuilder and null in python. I want to include Merge tool in my model. However:

at run 1, I want to Merge shapefiles s1 and s2
at run 2, I want to
Merge shapefiles s1, s2 and s3
at run 3, I want to Merge shapefiles s1, s2, s3 and s4 etc.

How can I include this to my modelbuilder model? I know how to make it manually, but how to make it in the model?

I can't just set the "Merge" as Parameter, and I don't see to option to make it via "Make Variable".

How the model looks like:

Exported python snippet :

arcpy.InfestationGradient("C:/Default.gdb/bb05","C:/Default.gdb/bb06","C:/Default.gdb/forest05","C:/infestGradient.gdb/buff_bb_t2","C:/infestGradient.gdb/buff_forest_t2","C://Default.gdb/bb05_Merge")


texts:
How do I check spatial reference on all parameters during Tool Validation?  Just looking for recommendations.  Here is what I have already. This must work with existing validation.

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()
    self.fcfield = (None, None)

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parmater
    has been changed."""
    if self.params[0].value and self.params[1].value:
      fc, col = str(self.params[0].value), str(self.params[1].value)
      if self.fcfield != (fc, col):
        self.fcfield = (fc, col)

        # Get the unique values of the field 'col' in the feature class 'fc'
        unique_values = [str(val) for val in
                            sorted(
                                set(
                                    row.getValue(col)
                                        for row in arcpy.SearchCursor(fc, None, None, col)
                                )
                            )
                        ]

        # Assign the unique_values list to parameters 2 and 3
        self.params[2].filter.list = unique_values
        self.params[3].filter.list = unique_values

    # Set the default values of parameters 2 and 3 to the first item in the list
    if self.params[2].value not in self.params[2].filter.list:
      self.params[2].value = self.params[2].filter.list[0]
    if self.params[3].value not in self.params[3].filter.list:
      if self.params[3].filter.list[0] != self.params[2].value:
        self.params[3].value = self.params[3].filter.list[0]
      else:
        self.params[3].value = self.params[3].filter.list[.1]

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    if self.params[2].value == self.params[3].value:
        self.params[3].setErrorMessage("Values cannot be the same")


texts: Representing data on thematic map in only two colors in qgis

texts:
I am very new to qgis, and i am doing a project on representing child immunization on thematic maps.

I want to represent the data which is even number in one color and odd number data in another color, i tried representing it in a thematic map.

In the style tab of layer properties I selected it as categorized and then selected the column which I want to represent it on map and selected color in color ramp and then classified it,but I don't know how to show even data in one color and odd data in one color.


texts: How to get FeatureClass or FeatureLayer from FeatureClassID?

texts:
I did a network analysis with ITraceFlowSolverGEN and the method FindFlowElements(). Then I did a queryId's with INetElements which returns the FeatureClassId and OID.
Now I want to get the FeatureClass Or FeatureLayer or anything that I can get the shape of the object.


texts:
I am calculating the VCI for Europe, but the ROI is too large for a single run in GEE. How can I divide it and process it in chunks so that when downloaded they fit into a single image/composite?
var bbox = ee.Geometry.Polygon([
  [-12.8,35.1],
  [-12.8, 71.5],
  [49.5,71.5],
  [49.5, 35.1],
  [-12.8,35.1]]);

Map.centerObject(bbox);
Map.addLayer(bbox, {color: 'FF0000'}, 'geodesic polygon');

// Assign a common name to the sensor-specific bands.
var LC8_BANDS = ['B2',   'B3',    'B4',  'B5',  'B6',    'B7',    'B10']; //Landsat 8
var LC5_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7',    'B6']; //Llandsat 5
var STD_NAMES = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'temp'];

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterDate('1990-1-1', '2000-12-31')
                                                      .filterBounds(bbox);// Landsat 8
var l5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA').filterDate('1990-1-1', '2000-12-31')
                                                      .filterBounds(bbox); //Landsat 5
var withCloudiness = l8.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: bbox,
    bestEffort: true,
    scale: 30,
  });

  return image.set(cloudiness);
});

var filteredCollectionl8 = withCloudiness.filter(ee.Filter.lt('cloud', 20)).select(LC8_BANDS, STD_NAMES);

var withCloudiness = l5.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: bbox,
    bestEffort: true,
    scale: 30,
  });

  return image.set(cloudiness);
});

var filteredCollectionl5 = withCloudiness.filter(ee.Filter.lt('cloud', 20)).select(LC5_BANDS, STD_NAMES);

var lall = ee.ImageCollection(filteredCollectionl5.merge(filteredCollectionl8));

var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['nir', 'red']).rename('NDVI');
  return image.addBands(ndvi).toFloat();
};

var ndviColl = lall.map(addNDVI);

// Define a combined reducer for min and max - one pass through the data.
var minMaxReducer = ee.Reducer.min().combine({
  reducer2: ee.Reducer.max(),
  sharedInputs: true
});

// Calculate regional min and max and add as image properties.
ndviColl = ndviColl.map(function(img) {
  var minMax =  img.reduceRegion(
    {reducer: minMaxReducer, geometry:bbox, bestEffort:true, scale: 30});
  return img.set({
    roi_min: minMax.get('NDVI_min'),
    roi_max: minMax.get('NDVI_max')
  });
});

// Filter out min and max region reductions that are null.
var filtered = ndviColl
  .filterMetadata('roi_min', 'not_equals', null)
  .filterMetadata('roi_max', 'not_equals', null);

// Create VCI index.
var vci = filtered.map(function(img){
 var id = img.id();
 var max = img.getNumber('roi_max');
 var Index = img.expression("(NDVI-min)/(max-min)", {
    "NDVI": img.select('NDVI'),
    "min" : img.getNumber('roi_min'),
    "max" : img.getNumber('roi_max'),
 })
 .rename('VCI')
 .copyProperties(img,['system:time_start','system:time_end']);
 return Index;
});

vci = vci.filterDate('2000-1-1', '2000-12-31');

Export.image.toDrive({
  image: vci.mean(),
  description: 'imageToDriveExample',
  scale: 30,
  region: bbox,
  maxPixels: 900000000
});

// Display the result.
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
var vciParams = {min: 0, max:1}
Map.addLayer(vci.mean(), vciParams, 'VCI image');


texts: Is there open source GPS post-processing software?

texts:
Possible Duplicate:
  Post Processing GPS Data with open source software

Here at the office we have a couple of 15-year-old Trimble 4700 and 4800 units but now we don't have the software to do the post-processing. Trimble wants to sell a new license so we can continue making full use of this equipment.

Is there an open source alternative to GPS post-processing?

I found GPSTk and RTKLIB and I wanted to know whether someone has had some experience using these or other similar software.


texts: What is the correct cartographic projection to calculate a Voronoi tesselation?

texts:
I have a set of points with known geographic coordinates. I would like to compute a Voronoi tesselation from these points in order to obtain areas of influence. The geometrical area of each Voronoi cell is then to be used to weight measurements obtained at each point.

Albeit area being my main interest, the computation of a Voronoi tesselation is a distance based algorithm. Considering that the set of points of interest cover about 5% of the Earth's surface, distances calculated in a, say, equiareal projection will incur in a relevant error. Conversely, if the option is to compute the Voronoi tesselation with an equidistant projection, the resulting areas will in their turn have relevant imprecisions.

The correct thing to do seems to be: calculate the tesselation with an equidistant projection, send the resulting cells back to the geographic domain, then re-project them with an equiareal projection and finally compute the areas. But this seems all too convoluted, would there be an alternative process of calculating such a Voronoi tesselation with precision?


texts: Re-projecting ground control points of similar images (Georeferencing in QGIS)

texts:
I have several maps to digitize. The base maps were Canadian NTS map sheets, all for more or less the same region. The useful information has been traced on top with handwritten markings. The maps were large so they had to be photographed rather than scanned.

It seems to me that I should be able to re-use ground control points to some extent after georeferencing the first map. Certain internal features on each map should be the same relative distance from each other.

What I want to do is manually georeference map A, and then warp map A's ground points onto map B and correct their positions. It seems to me that this would save time compared to manually entering all new control points - clearly there will be some error in the transform, but if control points are near enough to their true positions I can just shift them into place instead re-entering them.

This should because the two maps are in the same projection. The images of the base maps are same subject to translation, scaling, rotation, and some unavoidable error.

So this is a lot like re-projection, I just want to reproject a set of control points instead of using control points to re-project a raster image.

I'm using QGIS. Is there a way to accomplish this in QGIS?


texts: Options for displaying PostGIS vectors in OpenLayers

texts:
I am very new to web-mapping, so this question may be a little misguided.

I want to have a very simple web interface using OpenLayers (or even Google Maps if that will be easier) which can display a number of point and polygon tables that are currently in PostGIS.

I want this to be as simple (and as low maintainence) as possible:

I don't want to have to create tiles (unless there is a very good reason), or have a massive, complex layer of middleware.
Performance doesn't need to be phenomenal, as this is for internal use and thus doesn't need to scale beyond a couple of concurrent users.
Styling of the vector layers is of minimal importance.
Open source tools are much preferred as I have almost no budget for this.

What is the recommended approach for this?


texts: Rule based styling in Qgis using a combination of attribute(get_feature()) and ILIKE

texts:
I Have 2 Layers:

I want to Filter Layer 2 by an attribute from Layer1 for some rule based styles. I'm currently using following expression to realize that:

"C1"  = attribute(get_feature(‘Layer 1', 'ID', '1'), 'C') or "C2"  = attribute(get_feature(‘Layer 1', 'ID', '1'), 'C') or "C3"  = attribute(get_feature(‘Layer 1', 'ID', '1'), 'C')… etc.

It does its job, but due to the amount of features (a few thousands), performance is quite bad.

So I made a new field called "CG" in Layer 2 using following expression in field calculator:

"C1"  + ' '+ "C2"  +' '+  "C3"  +' '+ "C4"   +' '+ "C5"… etc.

The Field “CG” I can filter using eg. “CG” ILIKE “%C1 %”

Question:
How can I now combine the attribute(get_feature()) and ILIKE expression to Filter the Field “CG” in Layer2 by an attribute of Layer 1 in Field “C” ?  An example would be awesome.


texts: Attribute table blank after exporting from successful Display XY

texts:
I loaded a .csv file into ArcMap 10.7. There are no spaces or weird characters in the field names. The fieldnames are less than 10 characters long. There are no blank fields. I am able to successfully Display XY Data (i.e. ArcMap recognizes the decimal degrees for Latitude and Longitude) and successfully create an Events layer. However, when I go to export the data (i.e., the points I have displayed on the map and can see), the resulting shapefile and attribute table is empty. I have tried selecting all the records and then exporting and I have tried creating a feature class from XY data, to no avail. I have never had this problem before.
Does anyone have any ideas what is going on?


texts: Creating polygon from current display extent using ArcPy?

texts: I am trying to create a polygon from current display extent using Python. I would like to use the polygon to clip features in the current display extent or set a boundary limit for similar procedures.


texts:
How can I split a road network into different links splitting at intersections??
The whole network is in one single layer.


texts: Fail to upload gpx to Garmin GPS (GPSmap 60CSx) from Linux using gpsbabel CLI

texts:
I am trying to upload some gpx data (GPX XML format) to a Garmin GPSmap 60CSx unit using the gpsbabel (v.1.5.4) command line interface for Linux (Ubuntu 18.04).

I followed instructions from the gpsbabel website https://www.gpsbabel.org/htmldoc-development/fmt_garmin.html. Specifically, to test the connectivity I tried to download data from the GPS device into a gpx file, using the command

gpsbabel -i garmin -f usb: -o gpx -F blah.gpx

This worked fine so I then tried to do the opposite using the command

gpsbabel -D1 -i gpx -f Olympos.gpx -o garmin -F usb:

but no tracks and/or waypoints appeared in the unit after the transfer was complete. The USB default unit is identified correctly

gpsbabel -o garmin -F usb:-1

returned

0 3352480812 292 GPSMap60CSX Software Version 4.00

then using usb:0 or usb: as output device and debug output enabled in command line I see no error :

GPSBabel Version: 1.5.4
options: module/option=value: gpx/snlen="32" (=default)
Unit:   GPSMap60CSX Software Version 4.00
ID: 292
Version:    4.00
Capability A10:
Capability A100: D110
Capability A201: D202 D110 D210
Capability A301: D312 D302
Capability A400: D110
Capability A500: D501
Capability A600: D600
Capability A601: D601
Capability A700: D700
Capability A800: D800
Capability A801: D801
Capability A900:
Capability A902:
Capability A903:
Capability A904:
Capability A905: D900
Capability A907: D907 D908 D909 D910
Capability A908: D911
Capability A914:
Capability A916:
Capability A917: D917
Capability A918: D918
Capability A1010:
Capability A1011:
Link_type 1  Device_command 0
Waypoint: Transfer 100 Type 110
Route:    Transfer 201 Header 202 Type 110
Track:    Transfer 301 Type 302
Waypoint type: 110
Chosen waypoint length 14
cet_util: Converting from "UTF-8" to "US-ASCII", done.

what could be wrong ?
Is transferring to the device supported ?


texts: Is there a quick way to clear all attributes from a layer but leave the polygons in place?

texts:
I have a Shapefile layer in QGIS 2.6 with several polygons, each of which has data in over 100 fields. I need to create a new layer with all the same polygons, but with all their data fields blank (set to 0, Null or empty depending on the field type).
Is there a quicker way of doing this than pressing delete on each field of each polygon one at a time, I would have to do this more than 1000 times this way.


texts: Shapefile error in QGIS, rivers of a map

texts:
I am a beginner with QGIS. This file was sent to me along with others by my professor. One of them doesn't work. I'm mapping an Australian river catchment (Daly), one of the shapefiles is for the territory and the other one for the rivers contained in the territory. This last one isn't recognized by QGIS.

The message:

  Invalid Data Source: D:\Desktop\Tesis\Para Qgis\Rivers.zip is not a valid or recognized data source.

The files:


texts: Seeking low-resolution DEM of California?

texts: I'm trying to create a 3D model of California (which I'll manipulate in Blender) at a fairly low resolution — just enough for large features (like Mt. Tamalpais, Coronado, etc.) to be legible.  Wherever I look, I find DEM's at much higher-resolution that would require a lot of stitching together and would likely tax my computer's capabilities.  Any ideas?  A bonus would be a DEM including Baja California too, but I realize most US sources won't have that data.


texts: Why computing the distance between two ellipsoidal coordinates needs to be iterative?

texts:
I have found internal Android code that calculates distance between two spheroidal coordinates. It is complex and it is iterative. Why? Can't this calculation be done in "single shot"?

The Java code below:

private static void computeDistanceAndBearing(double lat1, double lon1,
    double lat2, double lon2, float[] results) {
    // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
    // using the "Inverse Formula" (section 4)

    int MAXITERS = 20;
    // Convert lat/long to radians
    lat1 *= Math.PI / 180.0;
    lat2 *= Math.PI / 180.0;
    lon1 *= Math.PI / 180.0;
    lon2 *= Math.PI / 180.0;

    double a = 6378137.0; // WGS84 major axis
    double b = 6356752.3142; // WGS84 semi-major axis
    double f = (a - b) / a;
    double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);

    double L = lon2 - lon1;
    double A = 0.0;
    double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
    double U2 = Math.atan((1.0 - f) * Math.tan(lat2));

    double cosU1 = Math.cos(U1);
    double cosU2 = Math.cos(U2);
    double sinU1 = Math.sin(U1);
    double sinU2 = Math.sin(U2);
    double cosU1cosU2 = cosU1 * cosU2;
    double sinU1sinU2 = sinU1 * sinU2;

    double sigma = 0.0;
    double deltaSigma = 0.0;
    double cosSqAlpha = 0.0;
    double cos2SM = 0.0;
    double cosSigma = 0.0;
    double sinSigma = 0.0;
    double cosLambda = 0.0;
    double sinLambda = 0.0;

    double lambda = L; // initial guess
    for (int iter = 0; iter < MAXITERS; iter++) {
        double lambdaOrig = lambda;
        cosLambda = Math.cos(lambda);
        sinLambda = Math.sin(lambda);
        double t1 = cosU2 * sinLambda;
        double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
        double sinSqSigma = t1 * t1 + t2 * t2; // (14)
        sinSigma = Math.sqrt(sinSqSigma);
        cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
        sigma = Math.atan2(sinSigma, cosSigma); // (16)
        double sinAlpha = (sinSigma == 0) ? 0.0 :
            cosU1cosU2 * sinLambda / sinSigma; // (17)
        cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
        cos2SM = (cosSqAlpha == 0) ? 0.0 :
            cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha; // (18)

        double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
        A = 1 + (uSquared / 16384.0) * // (3)
            (4096.0 + uSquared *
             (-768 + uSquared * (320.0 - 175.0 * uSquared)));
        double B = (uSquared / 1024.0) * // (4)
            (256.0 + uSquared *
             (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
        double C = (f / 16.0) *
            cosSqAlpha *
            (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
        double cos2SMSq = cos2SM * cos2SM;
        deltaSigma = B * sinSigma * // (6)
            (cos2SM + (B / 4.0) *
             (cosSigma * (-1.0 + 2.0 * cos2SMSq) -
              (B / 6.0) * cos2SM *
              (-3.0 + 4.0 * sinSigma * sinSigma) *
              (-3.0 + 4.0 * cos2SMSq)));

        lambda = L +
            (1.0 - C) * f * sinAlpha *
            (sigma + C * sinSigma *
             (cos2SM + C * cosSigma *
              (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)

        double delta = (lambda - lambdaOrig) / lambda;
        if (Math.abs(delta) < 1.0e-12) {
            break;
        }
    }

    float distance = (float) (b * A * (sigma - deltaSigma));
    results[0] = distance;
    if (results.length > 1) {
        float initialBearing = (float) Math.atan2(cosU2 * sinLambda,
            cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
        initialBearing *= 180.0 / Math.PI;
        results[1] = initialBearing;
        if (results.length > 2) {
            float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
                -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
            finalBearing *= 180.0 / Math.PI;
            results[2] = finalBearing;
        }
    }
}


texts: Why can PostGIS create LINESTRING from the same values of coordinates of the POINT type?

texts:
So, I needed to create linear objects from point objects.

I split the objects into points with the following query:
CREATE TABLE roads_arc_var3 AS SELECT (ST_DumpPoints(geom)).geom FROM roads_arc_var1;
I created the following table, so that it created two identical fields with the same values of coordinates:

CREATE TABLE roads_arc_var4 AS SELECT a.geom AS a, b.geom AS b FROM roads_arc_var3 AS a, roads_arc_var3 AS b WHERE a.geom = b.geom

I added the field "geoline" to the table:

SELECT AddGeometryColumn('roads_arc_var4', 'geoline', 4326, 'LINESTRING', 2);

I changed the values of the initial and final coordinates in some records in order to create lines from them, and I expected that in the 'geoline' field line values will appear only in those records in which the values of the coordinates in the starting point are not equal to the values of the coordinates end point, the rest will remain empty

UPDATE roads_arc_var4 SET geoline = ST_MakeLine(a, b);

But, the picture was different, see figure

What is the explanation for the PostGIS developers of this situation, because it's about 2D?

Or it is a problem with rounding off the values of coordinates, but in this case the user needs to understand what is its length, etc.


texts: Automatic extraction of manually selected points in ArcMap

texts:
I have an Orthomosaic and a DEM (both are in Geotiff formats and have RGB values) of a site, and I want to check their accuracy. For this purpose, I measured the 3D location (with mm accuracy) of some fixed features on-site with a total station. For accuracy comparison, I need the 3D coordinates of those fixed features from the orthomosaic (X,Y) and the DEM (Z).

I'm using ArcMap 10.1, and so far, the only way I could extract the above mentioned information is by placing the cursor on each feature and manually read the (X,Y) coordinates from status bar.

I have over 400 features to compare and manual reading of the coordinates takes a lot of time.

1) Is there a way that I could use click and save (or similar automated procedure) to obtain the (X, Y) coordinates from orthomosaic.

2) Also, I want to extract the Z coordinates from the DEM for a set of (X,Y) values. Is there an automated sampling procedure for this in ArcMap?

Note: I have the licenses for most 3D tools in ArcGIS including 3D analyst.


texts: How to filter out layers with QgsMapLayerComboBox?

texts:
I am using Qt Designer 4.8.5 to develop a plugin (with QGIS 2.14.3-Essen).

I recently started to use QgsMapLayerComboBox widgets instead of the conventional Combo Box widget. However, I want to filter out a number of layers so that only a few are shown.

When I used the Combo box widget, here is the code:

root = QgsProject.instance().layerTreeRoot()
main_group = root.findGroup('Main group')
layers = []
for group in main_group.children():
    for layer in group.children():
        layers.append(layer.layerName())
self.dockwidget.combo_box.clear()
self.dockwidget.combo_box.addItems(sorted(layers))

How can I translate the above code to work with QgsMapLayerComboBox? I have seen the QgsMapLayerComboBox.setExceptedLayerList() but not exactly sure how to create such a list.


texts: Applying same change to multiple map documents?

texts:
I have several map documents (around 30) that reference all the same files and display them all the same way. I created them by setting one map how I wanted it, then using "save as" to create a new map then adjusting the data frame extents in the layout. Of course now that I have them all created, I need to make the exact same change to all of them. Is there a more efficient process than just manually repeating the same steps in each document? If not, is there a better way I could have set up the project to facilitate changes like this?

I'm using ArcMap 10.0.


texts: MapBox GL geoJSON markers disappear when custom style loaded

texts:
I am using this example to load geojson points to to a mapbox basemap.

I view it on a browser locally and it works.

When I change the style reference to my custom style, the style will display, but the geojson entries will not. I think its trying to access a sprite(?) reference, but i'm not sure how to set that up. Style reference code:

var map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/elgirafe/cijvtamnh006uqtkqwqevcvxv',
    center: [-77.04, 38.907],
    zoom: 11.15
});

I tried to download my style, an change the sprite reference to "mapbox://sprites/mapbox/bright-v8" but that did not work either.

Can anyone shed some light on this?


texts:
Using vb.net, my custom tool needs to create a data table in my geodatabase to store some information for later use.  The user will be in ArcMap 9.3.1 with the geodatabase up and layers displayed.  I want to create the table at the same level as the dataset.  I can accomplish this with the following:

    gp = CreateObject("esriGeoprocessing.GPDispatch.1")
    gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    gp.createtable_management("C:\Dan\Data\geodatabase.mdb", TableName, "", "")

However, the problem is I cannot figure out how to have the geodatabase be something other than hardcoded.  How do I get the current geodatabase name and use a variable populated with that name in the createtable_management line?

I've searched high and low for two days now with no luck.

Thanks,

Dan


texts:
I've got this piece of code below that is supposed to display two base layers but now the OSM_LAYER disables(becomes blurry and inactive in the layer tree) the WORLD_LAYER, I'm not sure what the problem is. If it fits, I'd welcome the advice on how to change the OSM_LAYER to be in the format of WORLD_LAYER.

Heron.options.map.layers = [

 new OpenLayers.Layer.OSM(),

new OpenLayers.Layer.WMS(
        "World image",
        'http://www2.demis.nl/wms/wms.ashx?WMS=BlueMarble',
        {layers: "Earth Image", format: 'image/png'},
        {singleTile: true,
        isBaseLayer: true,
        visibility: false,
        noLegend: true}
)]


texts: Is there a way to speed up gIntersections for large spatial objects?

texts:
I'm attempting to perform a gIntersection of a big grid on a shapefile for country boundaries, so that I can calculate the % area of each grid cell that lies within the country. It takes a few seconds for smaller countries, but takes ages for larger ones.

I have a method that cut down the time substantially, but it still takes an unreasonable amount of time for larger countries. Here's a reproducible example, that uses gIntersects to subset the grid, which cut down on time substantially.

#### demo script for stackoverflow ####
# set up demo grid
lat <- seq(from = -11, to = 55, by = 0.5)
long <- seq(from = 70, to = 150, by = 2/3)
demo.grid <- SpatialPoints(coords = cbind(rep(long, times = length(lat)), rep(lat, each = length(long))))
demo.grid <- SpatialPixels(demo.grid, tolerance = 0.000001, proj4string = CRS("+proj=longlat +datum=WGS84"))
demo.grid.poly <- as.SpatialPolygons.GridTopology(getGridTopology(demo.grid), proj4string = CRS("+proj=longlat +datum=WGS84"))

# looks like a grid
plot(demo.grid.poly)

# get my file - I used the SpatialPolygonsDataFrame for level 0 administrative boundaries for China, downloaded from www.gadm.org
CHN <- get(load(file.choose()))
plot(CHN, add = TRUE, border = "blue") # China is located in the right place

CHN.p <- spTransform(CHN, CRS("+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
demo.grid.poly.p <- spTransform(demo.grid.poly, CRS("+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

plot(demo.grid.poly.p)
plot(CHN.p, add = TRUE, border = "blue")

# looks good

system.time(subset <- gIntersects(CHN.p, demo.grid.poly.p, byid = TRUE)) # find the grid cells that intersect with China
summary(subset)
demo.grid.poly.subset <- demo.grid.poly.p[subset] # make new set of gridcells that are subsetted to only those that intersect with China

plot(demo.grid.poly.subset)
plot(CHN.p, add = TRUE, border = "blue")
# still looks good
system.time(intersection <- gIntersection(CHN, demo.grid.poly.subset, byid = TRUE)) # create intersection SPDF. Go get sandwich, sleep, exercise, take a small vacation. Learn piano.

This is related to a question I asked a while ago (https://stackoverflow.com/questions/16918767/gintersection-on-very-large-spatial-objects). I'm facing a similar issue with a new dataset, and have a new method for solving the problem that has a new problem.


texts:
In an attempt to create an architectural site model from GIS data, I am trying to export a project from QGIS to .dxf, with the intention of setting up a cut file in Rhino.

QGIS tells me the export is successful, but when I open the exported file in Rhino, I can see all of the layers in the original QGIS file represented, but there is no vector information in the drawing.

Dxf files that I export are all about the same size, no matter how many layers I select in the export dialog, so I assume I am not checking a box somewhere. . .

I can get a single layer at a time out of this file using the "save as" function in QGIS, but when I try to coordinate them in Rhino, they appear to be based on different coordinate systems--importing at different sizes, and in different relationships to 0,0,0.

I have tried to check that they are all on the same CRS, and they appear correctly in QGIS, but they come into Rhino in completely different conditions. I have also made sure that all of the raster/shp files in the project have been converted to vector in QGIS, and again, this all appears correct in QGIS.


texts: Batch convert dgn files to separate shapefiles in ArcGIS

texts:
Actually, currently we are doing dgn -> shp file conversion using ogr2ogr. Batch convert dgn files to shp files using ogr2ogr

The problem is, at times the point layer is not getting read properly & few attributes are missing. But, if we convert the dgn using xxx.dgn -> Export -> To Shapefile (multiple) all the attributes of the features are getting retained and annotaion layer & point layer is getting saved separately.

But we are doing it file by file. Is it possible to run this as a batch using Arcpy ?

I've tried with the code given in this site https://geonet.esri.com/thread/56828

It work's perfectly. But the files are saved in a gdb. I want it to be saved as separate shapefiles.

import arcpy
import glob
import os
gdb = "D:/Arc/new.gdb"
arcpy.env.workspace = gdb
arcpy.CreateFileGDB_management("D:/Arc", "new.gdb")
reference_scale = "1500"
for file in glob.glob(r"D:\Arc\*.dgn"):
outDS = arcpy.ValidateTableName(os.path.splitext("d" + os.path.basename(file)) [0])
arcpy.CADToGeodatabase_conversion(file, gdb, outDS, reference_scale)


texts: How kriging variance is calculated in R gstat?

texts:
How does exactly the krige function (a wrapper to gstat and predict functions) from package gstat calculate kriging variance (minimized estimation variance) in Ordinary Kriging?

I wanted to reproduce in R results from How to make a prediction in Kriging using a semivariogram? (prediction equal to 6.88 and kriging variance equal to 3.14). Here is the R code for the same set of points and covariance function.

library(gstat)
library(sp)

p = SpatialPoints(cbind(c(0,0,3),c(0,2,1)))
p$z = c(10,7,3)
krige(z~1, p, SpatialPoints(cbind(1,1)), vgm(psill=3, nugget=1, range=6, "Exp"))

#Results
[using ordinary kriging]
  coordinates  var1.pred var1.var
1      (1, 1)  6.813431  2.002507

The prediction result is close (6.88 versus 6.81), but the kriging variance is very different (3.14 versus 2, even though the variance unit is squared which partially explains a greater difference).

The example I used in the linked post was made up, yet I also have tested examples from two other sources and was not able to reproduce results with krige as well (though results were a little more proximate than mine). What am I missing?

The equation used to calculate the kriging variance in the example was:

σ²ε = sill - [w1 ... wn λ] [C10]
                           |...|
                           |Cn0|
                           [1  ]

where σ²ε is the kriging variance, sill is the variogram sill parameter, wn the kriging weight of sample point n , λ is the Lagrange multiplier, Cn0 is the covariance between sample point n and prediction point.

A bonus question is, should the predicted values (6.88 versus 6.81) have been the same as well?

Looking at krige and predict source codes did not help because they are over my head.


texts: ArcMap query layer with Oracle Locator displaying no results

texts:
I have an Oracle 11g database setup with Oracle Locator.

I have created a table with the this query:

CREATE TABLE geotest (id NUMBER(38) NOT NULL, shape SDO_GEOMETRY);

Then I inserted one point with the this query:

INSERT INTO geotest VALUES(1,SDO_GEOMETRY(2001,3395,SDO_POINT_TYPE(-111.949439,40.722283,NULL),NULL,NULL));

The result is no data in the query layer in ArcMap. It allows me to create the layer, and validates. But then it makes me go through the advanced properties which I suspect means it can't see any of the rows to automatically pull the SRID.

Any guidance would be appreciated!

Update: So when I insert the metadata into the table using this:

insert into user_sdo_geom_metadata values ('geotest','shape',sdo_dim_array(sdo_dim_element('X',-180,180,0.005),sdo_dim_element('Y',-180,180,0.005)),4326);

And creating the spatial index, I can see the following point inserted into the database using this query:

insert into geotest values(1,'test',sdo_geometry(2001,4326,sdo_point_type(0,0,0),null,null));

So it would appear that a point with all 0s work, but anything beside that doesn't. I haven't found anything in the documentation that relates to this, or even any users complaining about that. Any either tricks that may be useful?

Update 2: Alright it appears that after I drop the index, and create it again, the points then show up. I am not an Oracle expert, so this confuses me, anyone have an explanation for this behavior? I thought once an index was created, as soon as a point is inserted, the index updates.


texts: Iterate through domain names and delete domains using ModelBuilder?

texts:
I have a list of unused/unassigned domains in table format. It's actually a table inside a GDB. I want to build a model that iterates through this table and deletes the domain.

Is this possible with ModelBuilder?

I tried to create a field value iterator, however, it won't let me link it up with the Delete Domain tool.

I know it's possible with python, but has anyone successfully attempted this with ModelBuilder.


texts: Feature selection list missing in editor Attribute window ArcGIS

texts:
Using ArcGIS desktop 10.7
The Attributes window in Editor identifies the selected feature. However, it is no longer showing the feature selection list. (The drop down menu on the right side of the feature field which allows for changing the feature record to available feature types from the attribute table).
How do I revert the editor settings to get this feature selection menu back? The list of feature types is also missing in the Create Features window. Appears editing features went missing when I added a different layer into my map document and now missing in any map document.


texts: batch georeference based on link file, output raster too big and in wrong place

texts:
I have a set of images that I need to georeference. All images have the exact same geographical coverage.
I have used the georeferenceing tool bar to create a link file.

I then used the script suggested in Batch Georeferencing in ArcMap/ArcPy?
to loop through the images.

arcpy.env.workspace = r"inputFolderPath"    # Folder with rasters to be georeferenced"
arcpy.env.outputCoordinateSystem = r"pathToLinkTable.prj"   #output projection
OutF = r"pathToOutputWorkspace"                         # Output folder
Ext = ".tif"                            # Format of the output rasters
LinkFile = r"pathToLinkFile.txt"           # Your link file
Transf = "POLYORDER0"
Resampl = "NEAREST"

# Process
ListRas = arcpy.ListRasters()
for ras in ListRas:
    basename = arcpy.Describe(ras).baseName
    outpath = os.path.join(OutF, basename + Ext)
    arcpy.WarpFromFile_management(ras, outpath, LinkFile, Transf, Resampl)
    print "Georeferenced {} successfully".format(basename)

All images are processes, however the output images/rasters are not in the expected geographic location.

I cannot see where the issue is. It must be in the source points/map points, but I am using the link table... any ideas what I should try next?

I am running the script in IDLE, as the python command line crashes when I try to run anything (unrelated issue but relevant information).


texts:
I can get informations on, for example, a stadium

[out:json];
(
  way[leisure=stadium]["name"~"Picot"](48.69241932797025,6.205483675003052,48.69873617918758,6.218605041503906);
);

out body;
>;
out skel qt;

But is there a way to get its entrance too (osm tag : entrance=main) ?

I think there is something better than the "around" command, but what... ?