texts: How can I do operations with points inside a buffer in GIS

texts:
I have a dataset with around 4k schools for with I have their location (longitude and latitude) and enrolment which I have uploaded to a GIS map. I want to do an operation with the competing schools in a 1km radious for each school. The operation is as follows:

I have identified that I need the following steps:

Draw a buffer of 1 km around each school
Identify all the schools j in that area for school i
Calculate the distance between the school i and all the schools j in the area
Run the CI formula for school i
Add the resulting number to the data table
Repeat for all schools

However, I have been unable to go past the second step. How can I do this in a way that covers all the schools?


texts: 'Undefined reference to' in standalone C++ QGIS application

texts:
I'm currently using Qt Creator to port a standalone PyQGIS application over to C++.
I'm using the OSGeo4W installation for importing the QGIS library.
However, I can't get my program to compile, as I'm getting the following errors when initializing the QgsApplication:

These errors are occuring when trying to do the following:

I've added these lines in the .pro file to import the QGIS library:

My IntelliSense correctly finds all of the files when trying to #include them, and it also autocompletes the QgsApplication functions such as setPrefixPath, initQgis or exitQgis. However, it won't compile because these functions are 'undefined references'?
I am a complete newcomer to C++ so I have no idea what the issue could be.


texts:
I am trying to create and insert data from ESRI shape files. This is one of the many features in the shape file, you can copy paste it and try it yourself:

DECLARE @g geography;
SET @g = geography::STPolyFromText('POLYGON ((28.397422886819093 70.5432251792964, 28.570915861286892 70.512498993109844, 28.540819573464251 70.4456402666085, 28.486718688957893 70.433538619104056, 28.353077926613448 70.512911092846991, 28.397422886819093 70.5432251792964))', 4326);
SELECT @g.ToString();

If I try to create a geography object from this on SQL Server 2008, it fails with the following error:

  Msg 6522, Level 16, State 1, Line 2
  A .NET Framework error occurred during execution of user-defined routine or aggregate "geography":
  System.ArgumentException: 24200: The specified input does not represent a valid geography instance.
  System.ArgumentException:
     at Microsoft.SqlServer.Types.SqlGeography.ConstructGeographyFromUserInput(GeoData g, Int32 srid)
     at Microsoft.SqlServer.Types.SqlGeography.GeographyFromText(OpenGisType type, SqlChars taggedText, Int32 srid)

What is wrong with this shape?


texts: How to leverage packages like Geocortex for small business GIS professionals?

texts:
Can anyone talk about how packages like Geocortex could be utilized by small business or working from home GIS professionals?

I would like to specifically know if Geocortex replaces any other packages or sets of tools that would make running your small business platform that much easier.


texts: Exporting vector to kml with ogr2ogr without extended data

texts: If I convert a shapefile as kml in QGIS, GRASS or ogr2ogr directly, defining the name and/or description field with the -dsco tag, all columns in the attribute table of the input layer are exported to the kml as ExtendedData. Is there a way to create the kml with the Name and Description only, i.e., without the Extended data?


texts: ArcGIS Online feature class error accessing from FME

texts:
I have a survey123 application. I am trying to connect the feature class connected with this application in ArcGIS Online to my local database. For this I am using FME Workbench. I have previously connected ArcGIS online feature class using FME Desktop and it worked fine. This time however I am getting the following error.
ArcGIS Online Feature Service Reader: Could not find a Feature Service item matching the Feature Service directive

I am using FME 2020.1
When I added the AGOL reader and chose this feature class I am getting column names in reader which means it is able to connect to feature class, when I press the see source data on reader I can see the data as well. It gives me error when I run the FME even if I don't have any transformer and just inspector.


texts:
I am using ArcMap 10.1 with a basic license.  It is my impression that the attachment feature has been upgraded to an advanced license.  My goal is to be able to add documents (ie. Land DEED) to a feature layer (parcel layer).

Is there any other work around that some knows of?


texts: Ogr2ogr: ERROR 1: PROJ: pj_obj_create: Cannot find proj.db

texts:
I start using GDAL 3.0.0 using conda build. I successfully installed GDAL on my laptop (Windows 10) and I am trying to ingest MapInfo tab file into Elasticsearch index. I understood that I have to create mapping first before ingesting it into Elasticsearch.

ogr2ogr -progress -lco MAPPING="D:\mymappings\map.text" -f "ElasticSearch" http://localhost:9700 "D:\MIProTab.tab"

But when I try to execute the above command then it is giving me the following error:

C:\Users\ABCD>ogr2ogr -progress -lco MAPPING="D:\mymappings\map.text" -f "ElasticSearch" http://localhost:9700 "D:\MIProTab.tab"
ERROR 1: PROJ: pj_obj_create: Cannot find proj.db
ERROR 1: PROJ: createGeodeticReferenceFrame: Cannot find proj.db
ERROR 1: PROJ: proj_as_wkt: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_wkt: Cannot find proj.db
ERROR 1: PROJ: createGeodeticReferenceFrame: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db
ERROR 1: PROJ: proj_as_wkt: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_wkt: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_wkt: Cannot find proj.db
ERROR 1: PROJ: proj_as_wkt: Cannot find proj.db
ERROR 1: PROJ: proj_as_wkt: Cannot find proj.db
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\Users\ABCD>

I checked other's solution where it is mentioned to have gdal.dll to be present on the same directory where ogr2ogr is present. Also, I have set the environment variable "PROJ_LIB" but all went in vain.
What is wrong here?


texts: How to refresh a layer in leaflet.js without flickering

texts:
We have an application with multiple layers, one of which is has to dynamically change based on new tiles from the server.
The base layer doesn't change, the overlay layer uses transparent png tiles, which are regularly updated.
The overlay layer has some code to make sure the tiles are not cached, without this the image doesn't change:

var dynLayer = L.tileLayer('base/dyn/{z}/{x}/{y}.png?{test}', {
    minZoom: 2,
    maxZoom: 5,
    tms: true,
    noWrap:true,
    test: function() {
        return Math.random();
    }
});

To refresh the changes the layer redraw method is called on a timer:

dynLayer.redraw();

Calling reload does not reload the images from the server.

This is based on an existing Windows application that's graphics heavy and has been modified to create tiles instead of drawing on the screen,

Using a vector layer isn't an option because it would require rewriting the server code.

The issue is that the overlay flickers when it is redrawn. Is there some plugin or other way of redrawing the layer without causing a flicker?


texts: Plotting of chart for BSI,UI,NDBI ,MNDWI,NDVI for an annual image under single chart in Google Earth Engine

texts:
I want to plot a chart for different index for an annual image under single chart only. But I do not know exactly what function to use.
I tried one but it is giving and error :
var bandsName = ['B01','B02', 'B03', 'B04', 'B05', 'B06', 'B07','B10','B11'];
var classNames = ee.List(['mndwi_2020', 'ndvi_2020', 'ui_2020', 'ndbi_2020','bsi_2020']);
// Define chart customization options.
var options = {
  lineWidth: 2,
  pointSize: 4,
  hAxis:{title: 'Bands'},
  vAxis: {title: 'Reflectance'},
  title: 'Spectral Profile for 2020 year indexes'
};

var spectraChart = ui.Chart.image.regions(
    image_2020, geometry, ee.Reducer.mean(), 30, classNames, bandsName)
        .setChartType('ScatterChart')
        .setOptions(options);
// Print the chart.
print(spectraChart);

Code link


texts: I am using ArcMap 10.2 and have a shapefile with parcel data.  I have all of the polygons labeled with attribute text.  How do i export the text to CAD?


texts: Using polygons made by Generate Service Areas as input feature class in ModelBuilder?

texts:
In ModelBuilder, I used "generate service areas" tool and made service area polygons.

I want to use them as input feature and connect to other geoproccessing tool, but I cannot.

Service area polygons remains as white　and I get:

  Error Message 000918

what should I do ?


texts: In ArcGIS 10.1 using Add-In, IApplication.Shutdown works fine on my PC, but when I install it on other computer it give error, something like "Not enough rights to close the application". The difference is: my PC - Win8, other - Win7 + in the domain. My best guess was that it is due to lack of Administrative rights that add in can't close ArcMap. Is there any workaround?


texts:
I have a database on PostgreSQL/PostGIS (9.5/2.3), and even if the final concept is to use it in QGIS (software that I have a good knowledge about), for now most of the people in the office are working with various versions of AutoCAD Map.

I defined a frame in which only the datas relevant to the area of interest are being loaded in QGIS through the use of views, and rules to allow their "updatability".

I would like to be able to do the same in AutoCAD Map, but it seems that when I load a view, it is automatically non-updatable, whereas with a regular table the update is possible. Does anyone know why it is like that, and, even better, how to do what I am trying to achieve?


texts: How to resolve case when routing result does not contain start and end point?

texts:
I have routing case shown on picture, SP=start point, EP=end point.

A-----SP-----B------------C----------D--------E-----EP------F

Routing result is BC, CD and DE. My points are not on route. I think that line direction is the key of problem.
If line AB has direction AB than my SP will be on route and source will be in point A, but if line AB has direction BA than my start point is not on route and source is in point B.

How to determine line direction to draw by myself segment SP-B?


texts: Configuring popup to automatically redirect to webpage using ArcGIS Pro / Online?

texts:
Is there a way to redirect the user to a specific webpage on clicking a certain polygon on a map?

Instead of a popup, the user must be sent to a website.

If this is possible how can this be achieved?


texts: Is there a way to route different vehicles within OSM2PO? I would like to start the service with a config file and then turn different flags on/off depending on whether the route is for a bike or a car or a pedestrian.


texts:
Problem

The picture is just an example of my case, but I am left with blank space on the left when I load all the tiles and I get error message in the console:

Failed to load resource: the server responded with a status of 404 (Not Found)
http://path/to/the/tile//base/8/6/53.png
etc.

I could guess and limit the extent of the layer, but I would really like to know is there any way to calculate it and also be able to get rid of the error message.

What I know about the tiles?

Originally they are served as TMS
They represent png images 256*256
They are in projection: EPSG:3857
In OpenLayers2 they were called like TMS layer
In OpenLayers3 (v.3.8.2) they are called like XYZ layer (that includes -y):

var baseLayer = new ol.layer.Tile({
  source: new ol.source.XYZ({
    url: 'http://path/to/the/tiles/base/{z}/{x}/{-y}.png'
  })
});

Thoughts and questions

How can I calculate extent of the tiles? Or in another words how can I find out zmin zmax xmin ymin xmax ymax or whole range of x/y?
What does x/y really represent? Are they coordinates of the left bottom corner or center of the tile?
In which direction does x/y increase?
How can I figure out x/y of at least one of the tile, for example left bottom one?

(I have modified the question for better understanding of the problem)


texts: Autocomplete not working for ArcPy in PyCharm

texts:
I'm trying to set up auto complete for ArcPy in PyCharm, but while submodules seem to have auto complete suggestions, the ArcPy module itself doesn't have completions provided.

Example of unsuccessful auto complete:

But when typing arcpy.gp. I see suggestions, as expected.

I thought that the submodules and functions/classes directly available from arcpy would be shown or help for specific functions would display.

I followed the instructions for setting up PyCharm to reference ArcPy. I ensured that the project interpreter was configured correctly, and also my PYTHONPATH in the application. I have also tried invalidating the cache and restarting PyCharm but did not see any improvement.

Is this auto complete supposed to work, and if so, how do I configure PyCharm to do it?

I am using PyCharm Community Edition 2016.2.3 and Arcpy as included with ArcGIS Desktop 10.3.


texts:
I know that the relationship of geoid height N and ellipsoid height h to orthometric height H is defined by H = h - N, but what if I have some position in Earth-Centered-Earth-Fixed (ECEF) coordinates?

Specifically, I am looking to relate the height of the geoid (for some particular position on the geoid) to the earth-centered height. (Using WGS84 here.)


texts: GDB lock remains when a layer is removed from the MXD

texts:
I have a script which removes a layer from the MXD, and then re-creates the feature from scratch in a GDB and adds the layer back in.

If I remove the layer manually before running the script, it runs fine. But if I remove the layer in the script, the script later fails when trying to re-create the GDB due to a lock. This is resolved by reloading the MXD.

  Traceback (most recent call last):   File "D:\RouteSmart_Work\Tools
  and Utilities\Python - Street Assignment
  Visualiser\StreetAssignment.py", line 60, in
      arcpy.CreateFileGDB_management(outputFolder, "StreetAssignments.gdb")   File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 16594, in
  CreateFileGDB
      raise e ExecuteError: ERROR 999999: Error executing function. Failed to execute (CreateFileGDB).

I've tried implementing arcpy.Compact_management (credit to StacyR from another and gc.collect() but these still don't seem to clear the lock. Could anyone advise if I'm missing anything? The script is meant to be run frequently, without the need to remove the layer each time.

def clearWSLocks(inputWS):
  '''Attempts to clear locks on a workspace, returns stupid message.'''
  if all([arcpy.Exists(inputWS), arcpy.Compact_management(inputWS), arcpy.Exists(inputWS)]):
    return 'Workspace (%s) clear to continue...' % inputWS
  else:
    return '!!!!!!!! ERROR WITH WORKSPACE %s !!!!!!!!' % inputWS

arcpy.env.overwriteOutput = True

servLoc_layer = arcpy.GetParameterAsText(0)
streets_layer = arcpy.GetParameterAsText(1)
outputFolder = arcpy.GetParameterAsText(2)

# Remove the layer if it's already in the data frame
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
if arcpy.Exists("FC"):  # Use name rather than layer object in case name doesn't exist
    lyr = arcpy.mapping.Layer("FC")  # Create layer object in order to remove it
    arcpy.mapping.RemoveLayer(df, lyr)
    arcpy.RefreshActiveView()
    del lyr, df, mxd

outGDB = os.path.join(outputFolder, "GDB.gdb") # Geodatabase

xyLineOutput = os.path.join(outGDB,"FC")
if arcpy.Exists(xyLineOutput):
    arcpy.Delete_management(xyLineOutput)

gc.collect()

if arcpy.Exists(outGDB):
    arcpy.AddMessage(clearWSLocks(outGDB))
    arcpy.Delete_management(outGDB)
arcpy.CreateFileGDB_management(outputFolder, "GDB.gdb")

EDIT:

In this instance I've instructed the user to remove the layer manually before running the script, to keep the one GDB as the workspace. It's only a small geoprocessing task, so the data management will be simple. I will bear Vince's comment in mind for the future.


texts: How to create polygon from bounding box coordinates (PostGIS/cartoDB)

texts:
I have data to import to cartoDB in CSV format in which there is a bounding box column holding two lat/long coordinates representing the corners of a bounding box, like

53.113206182344,-9.15629032156; 53.119089697211,-9.141677599209

I need to retain this information as text in its own column, but I also would like to derive a polygon and store its geometry as the_geom. I'm assuming that the polygon geometry would be created post import through a SQL update ... but need guidance in how to approach this.


texts: Are there any free KML to DAT(destinator compatible) converter available ? I tried to use global mapper but was in vain. Also used ITN converter but the line was converted to points.


texts: How can I modify the initial configuration of QGIS, using pyQgis?

texts:
I want to modify the initial configuration of QGIS.
I read what can be done with the environment variable PyQGIS_STARTUP, but I do not know how to do it.

What I want is to remove some QToolBar, but before I start QGIS.

I know I can remove it using:

main_window.removeDockWidget ([i for i in main_window.findChildren (QDockWidget) if i.windowTitle () == 'Panels'] [0])


texts: Select geographic transformation used when reprojecting

texts:
Is it possible to select the geographic transformation that is used when reprojecting the layers in QGIS?

It seems like ESRI defaults to EPSG 1900 when reprojecting the layer from NAD 83(HARN) to WGS 84 (which is weird considering EPSG 1901 supersedes it)

https://epsg.io/1900

Meanwhile, QGIS seems to default to EPSG 1580 when reprojecting the layer from NAD 83(HARN) to WGS 84 and does not give the option to change it.

https://epsg.io/1580

I believe this might have something to deal with selecting the geographic transformation in QGIS:

  +towgs84=0,0,0,0,0,0,0

And I believe that this stack exchange question is related to my issue:
Converting NAD 83 to WGS 84 in PostGIS


texts: Using processing using other output (result from another algorithm) as input

texts:
I have a problem. I am developing a standalone app in Python and I need to use processing algorithms.

I was able to use one algorithm without errors but I need to run other processing algorithm which input is the output of the first one. So, my program stops in the first algorithm. Why it doesn't read the second one?

CODE:

import sys

from qgis.core import QgsApplication, QgsVectorLayer
from PyQt4.QtGui import QApplication
from PyQt4.QtCore import QFileInfo
import GdalTools_utils
import ftools_utils
import osgeo
import ogr
import gdal
app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

from processing.core.Processing import Processing
Processing.initialize()

inputLayer = QgsVectorLayer(...)
inputMask = QgsVectorLayer(...)

extent_rect = inputMask.extent()

xmin = extent_rect.xMinimum()
xmax = extent_rect.xMaximum()
ymin = extent_rect.yMinimum()
ymax = extent_rect.yMaximum()
extent = str(xmin) + "," + str(xmax) + "," + str(ymin) + "," + str(ymax)

cellsize = 30
outrains10 = '...'

count = inputLayer.featureCount()

rain_10 = 'rains10'
days_10 = 'days10'

Processing.runAlgorithm("grass:v.surf.idw", None, inputLayer, count, 2.0,rain_10, False,extent, cellsize, -1.0, 0.0001, outrains10)

Processing.runAlgorithm("saga:clipgridwithpolygon", None, outrains10 , inputMask, outclip_rains10)

So the first algorithm grass is OK. The output is fine. But I want to use it as input in the second one, but the code stops there.

UPDATE

The code stops because I didn't know that I have to install SAGA to use in Linux. Now I have SAGA installed (I run from QGIS and everything is OK) but I still don't know how to configure SAGA to read through processing, such as the line:

Processing.runAlgorithm("saga:clipgridwithpolygon", None, outrains10 , inputMask, outclip_rains10)

UPDATE

I still can't get SAGA working but with many many tries, I found that, when I run the script, an error (Segmentation fault) and others are returned. So debugging my code I found that is the input raster which gives the error:

example:

grid = QgsRasterLayer('/.../raster.tif')

Does anyone knows why this is happen?? Can you help me please?


texts: Clipping 5 km area of watershed upstream from outlet using ArcMap

texts:
I am using ArcMap.
I have the delineated watersheds for more than 200 sites. What I need is to clip the 5-km area upstream of the sites. https://i.stack.imgur.com/S19Me.jpg.
I cannot do this manually for all sites and I am looking for any ideas that can be helpful to cope with the issue.  It would be doable even with rasters, knowing the size of raster pixels.
Is there any tool to select like 2000 pixels upstream a specific point?
It cannot be like a buffer as the buffer create a circle or square all around the sites not the upstream.


texts: Changing QGIS attribute table appearance

texts: How can I change the appearance of the attribute table in QGIS, so that I can view the cell's font an its size smaller in my computer's screen?


texts: Determine If PostGIS Geometry Is a Circle

texts:
I have a PostGIS database table that contains several generic geometries in one of it's columns. Sometimes I have geometries that are intended to be circles which I insert into the database as the result of a center point and a radius point using something like this:

INSERT INTO myTable (drawing)
VALUES (
    ST_Buffer(
        ST_MakePoint(0, 0)::geography,
        (SELECT
            ST_Distance(
                ST_GeographyFromText('POINT(0 0)'),
                ST_GeographyFromText('POINT(1 0)')
                )
            )
        )::geometry
    )

Now I would like to query the table for all the values inserted in this way, using something like this:

SELECT * FROM myTABLE WHERE FUNCTION_SAYS_GEOMETRY_IS_CIRCLE(geometry)

Does anybody have a recommendation? If this involves inserting my data differently that is okay. I have thought about adding an extra column for "is_circle" or something but that feels like a hack. I would prefer some insight from the geometry itself.


texts: I installed geoserver software on my system as well as the JRE(Java Running Environment), but the geoserver admin through which I am supposed to log in appears as 'short cut', and more over if I start the server, it runs for some minutes and closes. I installed different versions of JRE as well as geoserver software but the same thing kept occurring. Please is there something am not doing right or omitting?


texts: OpenTripPlanner API Elevation values unclear

texts:
I'm using the OpenTripPlanner  Planner API to find directions.

The API returns XML with the geometry of the line (2D) and the elevation data separately, for example

The 2D line geometry is as encodedPolylineBean in:

response
.plan
..intineraries
...legs
....legGeometry
.....points

And looks something like

wnseIhqzL??DCl@m@b@m@lA}AjB}B`@g@d@q@\\i@l@y@Vc@vAqCbAeCvAsDx@qBdCaG^{@DGr@_AxFiAjDcAv@]jAg@h@WlBcA\\Sz@c@@RPxA@`@Ad@Ex@@NDJBJ^dATp@h@xCBN^x@h@jAX^`AhApAzAbAfAfAhAtAxA|AfArBl@|Ad@vAXxCrAt@X^Vn@`@pA`ApA`AvAdAzAd@\\Ep@KFAbDgAhBwBnBkCBE`@i@p@_@t@a@zAa@VGzAw@x@jCn@pBFPNf@j@dBvAxDHTXz@lAfDvA|D?D`BfFz@hCFPz@nCdA~CQb@

But the elevation data is in:

response
.plan
..intineraries
...leg
....steps
.....elevation

And Looks like:

<elevation>
  <first>5.322</first>
  <second>73.87</second>
</elevation>
<elevation>
  <first>15.322</first>
  <second>74.25</second>
</elevation>
<elevation>
  <first>25.322</first>
  <second>74.56</second>
</elevation>
<elevation>

What I don't understand is:

What is the difference between the <first> and <second> elevation?
Why when the PolylineBean line is decoded does it have fewer points than the number of values given for the elevation?

For reference, a full XML response is available here


texts: How can I programmatically identify missing sidewalks in OpenStreetMap data?

texts:
I have a web map that is using the pedestrian routing from Mapbox's directions API. However, users have seen very long circuitous paths in certain locations because it appears that sidewalks (and/or foot access being allowed) are missing in the underlying OSM data. Is there an easy way to identify these gaps to fix within a major a city? Is there a way to query Overpass-Turbo to display roads/streets that have access:foot=not specified (or access, not food)?
Update
According to Mapbox

The Directions API (Walking profile) allows routing via all roads except:

motorways
trunk roads
all other roads where the value of either the 'foot' or the 'access' tag is set to any of 'no', 'private', 'agricultural', or 'forestry'

Roads other than motorways and trunk roads not having a 'footway', 'foot', or 'sidewalk' tag associated do not hinder routing as by default these routes are assumed to be accessible by pedestrians.


texts: Changing CRS of downloaded shapefile using QGIS?

texts:
I have looked through all the other questions which might have received answers to the question I am currently posting. However, after trying all their recommendations I am still stuck with the same issue.

I have downloaded data which provided shape files with are set within the coordinate reference system (CRS) as EPSG:4230 - ED50. I need this shapefile in the CRS of EPSG:4326 - WGS 84.

I have tried adding new vector layer with the downloaded data shapefile. I then saved the file with the CRS that I desire. When I upload it again it still is recorded as EPSG:4230 - ED50, does anyone know why I cant change it. I am using QGIS versions 3.2


texts: Visibility of multiple points from one other point?

texts:
I am fairly new to GIS and am struggling a little with one of my maps.

I have to see if 49 wind turbines (200m) are visible from a radar that is 39.6m high, I understand I need to do viewshed but have no idea what to put in each of the boxes. The area being used is a hilly area in Scotland so obviously this will have an impact.

I am using ArcMap 10.5 and so far have all of my turbines plotted and the radar.

What are the next steps to see if the turbines are visible from the radar?

I don't need the radar scope I essentially need to see if there is visibility between numerous points.


texts: Setting the study area boundaries in QGIS for nearest neighbor analysis?

texts:
I've made a cursory (i.e., two hour) search of the Stack Exchange and the Internet to find out how we can set the study area boundary for the nearest neighbor analysis command in QGIS 3.16. Sadly, I have found no definitive answers.
So here are my questions:

Is it possible in QGIS to set the study area boundaries? I have already tried running the NNA operation with my study area boundary on, my study area boundary off, then reprojected to WGS 84 (to make sure it met the lat longitude of the point data layer instead of relying on projections on the fly) and then by setting the scale of the window so that the study area just fit inside the window. In all four cases, I always got the same results. Sadly, there appears to be no option in the Nearest Neighbor Analysis window that let's us set the study area boundary.

If it is possible to set the study area boundaries for the NNA, how do we do it?


texts: Showing both sides of street from roads layer?

texts:
I am trying to show parking occupancy during peak hours and I need to offset a roads layer to show on both sides of the street.

How to get the lines to offset?

I am trying to make a map like the one shown below


texts: Extract row & column numbers of a raster at points in QGIS

texts:
Is there a tool in QGIS that will extract the row and column numbers of a raster at a set of points? All I have been able to find is the following:

Raster Analysis > Sample raster values
(Plugin) Point Sampling Tool
GRASS > Vector > v.sample

None of these return the row & column number of the point, as far as I have been able to tell.

I also found this question, but it is rather in-depth and I was hoping to find a tool which could do it automatically if possible.


texts: How to add a buffer in QGIS that updates itself

texts: I wish to add a buffer into project that will update the buffer fields when new points/ polygons are added into the project


texts: Comparison of open data portal solutions

texts:
For many organizations, such as local/county/state governments, providing public access to GIS data is a challenge.
Simply knowing what data is available, what it is, who owns and maintains it, and where to get it is probably the most immediate challenge. Therefore a GIS data inventory is a necessity. Of course, creating one is one thing, getting people to maintain it (and assume responsibility for their data) it is quite another. Therefore getting buy-in from internal stakeholders is also very important.
Another important consideration is a solution's support for metadata. If the data cannot be described, it is almost useless. With some solutions there may be limited support for metadata in the form of a limited number of proprietary fields. Support for metadata in its full and original form seems to be rare.
There may be limitations on the size, number of records/attributes, field lengths, or the geometric complexity of common GIS data, such as address points or parcel data that might make some solutions unusable. Often, the sheer size of some types of data (e.g. LiDAR, aerial imagery) makes it difficult or impossible to host them over the internet in their raw form, and instead they must be delivered via physical media (e.g. hard drives). At best, these large files may be hosted on something like an FTP server or Dropbox/Box, but it may be quite costly to do so.
Additionally, in many implementations, the spatial nature of the data may not considered to be important - or considered at all - making finding, viewing and using spatial data spatially difficult or impossible. Things like being able to symbolize the data in any useful way may be nonexistent.
Lastly, keeping the data up-to-date is vital. If there is no automated way to update data then the whole thing becomes unworkable for non-trivial amounts of data.
So with that out of the way, the main question is:

What open data portal solutions are available and how do they compare?

Side questions:

What are the most important considerations when deciding on an open data portal solution?

Are there scenarios where using more than one portal is beneficial (and cost-effective)?

Can spatial and non-spatial data be successfully integrated into the same portal?

Are there any examples of open data portal implementations you would consider to be successful or exemplary?

A list of some implementations is in List of US government free GIS data online as well as the Wikipedia article on geoportals.

Addendum:

A list of portal solutions that we are currently evaluating:

ArcGIS Online
ArcGIS Open Data
Portal for ArcGIS Server (on-premises)
Socrata and Mondara
Google Maps Engine and Gallery

Other portals we haven't yet considered:

DataHub
CKAN
DKAN
OpenGeoportal

Not really portals but work well for serving flat files:

Dropbox
Box
FTP (on-premises)

Rolling your own, either using hosted storage/servers or on-premises

In our case, because we have an ESRI ELA and have an ArcGIS-based spatial data infrastructure, it seems quite likely that at least one of the ArcGIS solutions will be selected.
Personally, I'd rather see us use more open source software, but it seems unlikely we'd be able to put something together that could compete with ESRI's offerings, at least not without a paradigm shift in the way people use GIS in our organization.
However, I hope that this question becomes a useful resource to anyone looking to improve access to public data.


texts:
I have performed Line Density and the results is a Raster...I want to make polygons from this Raster, but I cannot do it...

It shows me when I put this raster as Input raster:

"One or more Dropped items were invalid and will not be added to the control."

Does someone know how to do this job in ArcGIS 10.1?


texts:
Is it possible to connect ERDAS 2011 with tools of OpenOffice instead of Office tools? I need to send a map composition made in ERDAS to Presentation of OpenOffice. I don't have Word or Power Point and I don't want to install it. I use OpenOffice.

Or is there other any way to send/export my map composition from ERDAS to OpenOffice.

I sent my map to JPEG, it works but resolution is bad to print it out and even to have it as the image in jpeg.

Thanks a lot.


texts:
I'm making this loop, it's part of a script I'm making to illustrate a percentage in a polygon (let's say instead of a pie) but then I realized it looks good as an animation; when I start plotting, the polygon goes black every now and then, I want to be able to see the line between the two sections and with the black polygon I can't

library(units)
library(sf)
library(dplyr)
library(osmdata)
library(spatial.tools)

xalver <- osmdata::getbb("xalapa veracruz", format_out = "sf_polygon")
ea <- sp::bbox(as(xalver, "Spatial"))
ea2 <- ea
for(i in seq(.001, .015, by = 0.001)) {
  ea[4]  <-  ea[4] -i
  ea2[2] <- ea[4]
  oik <- bbox_to_SpatialPolygons(ea) %>% sf::st_as_sf() %>% st_set_crs(4326)
  oik2 <- bbox_to_SpatialPolygons(ea2) %>% sf::st_as_sf() %>% st_set_crs(4326)
  anka <- rbind(oik, oik2)
  st_intersection(xalver, anka) %>% plot()
#st_area() -> osh

#  if( drop_units(osh[1])/drop_units(osh[2]) < 0.5) break
  Sys.sleep(.3)
  }


texts: Line cracks in tiles generated with gdal2tiles

texts:
Objective: generate mercator tiles 256x256 from a huge jpeg2000 image.
Tool used: gdal2tiles.py
GDAL version: 3.2.2
Command line: gdal2tiles.py -z 19 -r [any] image.jp2 tiles/
Source image: image.jp2 is in CRS Lambert93 (2154), coming from IGN ORTHO-HR open database (resolution 20 cm per pixel) https://geoservices.ign.fr/documentation/diffusion/telechargement-donnees-libres.html#ortho-irc50cm-et-hr
Problem: there are some cracks in all tiles generated, visible on road or building edges (as shown in illustration). The artefact doesn't depend on the resampling method used.
Illustration:

The samples are 256x256 pixels tiles, generated with gdal2tiles, with different -r options. 1_original.png is an extract (done with QGIS) from jp2 original image, for the same area than the tile.
The source CRS is well defined in source image metadata, so I don't use -s option on command line to specify it (I tried, but it doesn't change anything as expected).


texts: Create points at changing line attributes

texts:
I have two layers:

Layer1: Linelayer with different types (stored als attribute values. Example: soil/material in which pipes are installed)

Layer2: Point-Layer to mark transition between different Line-types.

I wonder if it's possible to create the points automatically (one point at every transition of my attribute line-type):

Edit1:
A tool that just sets a point at the end of every line segment is not enough as there can be two or more linesegments with the same attribute in succession.
(See screenshot2: the point in the red rectangle should not exist in the point layer!)


texts: Removing small segments of geometry from shapefile using FME?

texts:
I am working on road dataset on a country level.

The gdb contains the entire road network which we have extracted from 3rd party source provider.

But the road network which we have fetched contains a lot of small dangling edges.

How can I remove all of such dangling nodes in bulk using FME?

I cannot do it manually, but it will take a lot of time and the sheer number of such dangling nodes are in lakhs.

Also I am not able to find the transformers like MRF2DCleaner and MRF2DDangleRemover as I am limited by the licencing level.


texts: ArcPy arcpy.gp.ExportMapServerCache invokes error: A secure connection with this site cannot be verified

texts:
What worked last week on many computers no longer works this week. My great ArcPy script gets a users current display extent and invokes gp.ExportMapServerCache to copy a web  cache locally. It was working very well across: 10.1 , 10.2 and 10.3 ArcGIS desktop.

When gp.ExportMapServerCache is reached in the script a dialog box is produced. The dialog boxes buttons are unresponsive. I can't get any other details about this apparent certificate information.

I don't have any issue adding the service to ArcMap and I can see it's details in the rest. This last week has seen a metric tonne of IT updates roll through. I guess somethings eaten something else.

How should I start diagnosing this?

Security Alert Dialog:


texts: The updateCursor doesn't update all the cells in the shp, just some of them

texts:
Don't know for sure what is happening, but when I use arcpy.UpdateCursor it only writes some of the lines and not all of the lines.

The code is:

cursor = arcpy.UpdateCursor("FInal1a1.shp")
for row in cursor:
    row.setValue("NomArq",n)
    cursor.updateRow(row)
    with arcpy.da.SearchCursor("Final1a1.shp",['Id','Diametro','NomArq']) as Buscador:
        with open(nombarch2,"a+")as c:
            for row in Buscador:
                for cell in row:
                    conv= str(cell)
                    c.write(conv)
                    c.write(" ")
                c.write("\n")
            f.close

And the output is something like this: (it's the third row)

2 75 21
2 75 0
3 100 0
3 100 0
3 200 0
2 75 21
2 75 21
3 100 0
3 100 0
3 200 0
2 75 21
2 75 21
3 100 21
3 100 0
3 200 0
2 75 21
2 75 21
3 100 21
3 100 21
3 200 0
2 75 21
2 75 21
3 100 21
3 100 21
3 200 21
2 90 32


texts: Query builder - IN function returning zero rows

texts:
I have a layer of several areas (polygons) each having a unique number, I want to make a new layer only containing some of the areas. I made a duplicate of the original layer and went to query builder and tried the following expression:
"n2000_nr" IN ('16','17','18','24','30','33','35','39','41','48','49','52','53','57','68','76','95','114','173','186','233')

But when I test it, it returns zero rows. I tried only selecting the 5 first areas and that went fine.
Is there a limit or what is the problem?


texts: GeoServer: Can't open Layer Preview for single layer

texts:
In Geoserver Layer Preview I can't open preview for single layers, in case of Layer Groups it works fine, but when I click single layer preview, pop-ups download window with a "wms" named file which includes the text below:
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://dev.ozorix.com/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1" >
<ServiceException code="LayerNotDefined" locator="layers">
Could not find layer mygeorgia:mg_5_roads
</ServiceException></ServiceExceptionReport>

I searched the issue in Google for hours but didn't succeed, tried to change styles, projection, edit wms links - no improvement.
Does anyone know how to solve the issue?


texts: How to reproject data in EPSG:4326 to EPSG:3857 using SharpMapv0.9 for display on Google Maps

texts: I have data stored in SQL Server 2012 in EPSG:4326. I need to display this data on a Google Map in EPSG:3857, but without an ST_Transform function, I cannot determine how to convert my data on the fly to EPSG:3857 to display them properly.


texts: I want to mosaic Landsat 8 scenes but the final product has dark areas on overlaps and on the background of each scene. Help to remove the dark areas.


texts: Testing quality of training data in Google Earth Engine

texts:
I need to manually define landcover training data in order to train a classifier for landcover change detection. My question is, if there is a way to test the correctness of the manual landcover classification. My issue is, that it is sometimes hard to detect differences in landcovers by eye (e.g. if a pixel is bare land, shifting cultivation or shrubland etc.).
Is there maybe an option to automatically multiply training points in order to test the classification manually afterwards?


texts: PostGIS populating existing table with imported shapefile

texts:
I have a PostGIS database with an empty table ch01.highways. With the following code, I want to populate the table:

INSERT INTO ch01.highways (gid, feature, name, state, geom)
SELECT gid, feature, name, state, ST_TRANSFORM(geom, 2163)
FROM ch01.highways_staging
WHERE feature LIKE 'Principal Highway%'

The ch01.highways_staging table was imported with the shp2pgsql-gui. For both tables, the owner is postgres.

When I run the command the ch01.highways remains empty, no data are imported.


texts: Generating map with slopes <10% using QGIS?

texts:
I am on QGIS 2.18.2, trying to generate a map that shows all areas with slope steeper than 10% in a specific color (red in this case) and hides the rest. I followed most of the instructions in this post: Slope Analysis in QGIS but I still have some problems, unfortunately.

Here is what I did so far: I got the raw data from the client in the form of points, imported them into QGIS and, per the instructions from the post above, generated a DEM layer using the interpolation tool. The result had a few weird lines in it. I used IDW as an interpolation method, because when I used TIN it looked even worse. Next I did a slope analysis, which resulted in this picture:

My QGIS is in German. The slope layer is called "Neigung". The weird-looking areas around the edges are expected, that data is missing from the raw file (but that's ok, the only data I need is the area of a municipality in the middle of the area). The problem is that I don't really know what unit the output is in. It should be degrees, but I highly doubt that the steepest slope in the area is a mere 4.56 degrees.

For the slope analysis, I left the z-factor at 1. The points and the raster layers are all in Gauss-Kruger, so everything is in meters, if I understood it correctly. And per the above post, the factor should be one if the vertical and horizontal units are identical, which they are in my understanding (the data points give height in meters and are spaced 25m apart in each direction.

I don't know what exactly the slope tool generated. It should be degrees, but the range does not really fit the situation.


texts:
I have been trying to use the azimuth expression function editor between 2 points.  However, I am not getting it to work.  I am obviously misunderstanding how to create the expression itself.  This is what I have

=  degrees( azimuth( make_point( "PrefRef"  'TR 15 NW 100' ), make_point( "PrefRef"  'TR 15 NW 114' )))

PrefRef being the point label and ID.

  Expression is invalid Error is  Parser Error: syntax error, unexpected
  EQ

  Eval Error: No root node! Parsing failed?

I am using QGIS 2.18


texts: Multiple prompts when logging in outside of domain

texts:
I work entirely with secured services in CMV and everything works fine when using them on a computer that is within my company’s domain.  Our IT has it setup so that you are prompted once to login to the IIS server then prompted a second time to login to AGS (not ideal/pretty but it works).

However, when I try and access the same web map on a computer outside my company’s domain I am prompted to login for every secured service that is within the operational layers of my web map.

Please let me know if anyone has info that will help me resolve the multiple prompt issue.


texts: WFS queries data through time field in GeoServer

texts:
Now using ol4.6.5, GeoServer uses 2.17.2, the following error appears

java.lang.RuntimeException: Parsing failed for Filter:
java.lang.RuntimeException: Unable to find function During Parsing
failed for Filter: java.lang.RuntimeException: Unable to find function
During Unable to find function During


texts: Python tool won't stop running-- Not an infinite loop issue

texts:
I'm having trouble with a python script.  It is essentially a near duplicate of another script which is running fine-- it is using different files, but apart from that the tools used and the order in which the tools are used is basically the same.

When I go run my duplicate tool, it won't stop running.  It doesn't have any logic in it which would cause it to go into an infinite loop.  It is essentially just a string of tools and outputs.  Does anyone know what it could be about this loop that is keeping it from successfully finishing?

I've tried including print statements, but it seems to freeze so the print statements won't print.

Here is the code (with path named edited):

import arcpy

arcpy.geoprocessing.env.overwriteOutput = True
#Variables
asOldParcelsFC  = r"string"
asBoundary = r"string"
scParcelsInAS = r"string"
newNotInOld_preclip = r"string"
newNotInOld = r"string"

arcpy.MakeFeatureLayer_management(asOldParcelsFC, "asOldParcels_layer")
arcpy.MakeFeatureLayer_management(asBoundary, "AltamonteBoundary")
arcpy.MakeFeatureLayer_management(scParcelsInAS, "scParcels")

arcpy.AddJoin_management("scParcels", "PARCEL", "asOldParcels_layer", "PARCEL", "KEEP_ALL")
arcpy.SelectLayerByAttribute_management("scParcels", "NEW_SELECTION", "asOldParcels.PARCEL IS NOT NULL")
arcpy.SelectLayerByAttribute_management("scParcels", "SWITCH_SELECTION")
arcpy.CopyFeatures_management("scParcels", newNotInOld_preclip)
arcpy.Clip_analysis(newNotInOld_preclip, "AltamonteBoundary", newNotInOld)

arcpy.Delete_management(newNotInOld_preclip)
arcpy.Delete_management("AltamonteBoundary")
arcpy.Delete_management("asOldParcels_layer")
arcpy.Delete_management("scParcels")


texts: Stretch image to fit boundary coordinates in TileMill

texts:
I would like to add an image as a layer to a Tile Mill map by specifying the coordinates. Is this possible? If so, how would I do this?

I've seen this being done, where an image is set to certain coordinates, but I would like to stretch an image to fit to two corner coordinates.


texts: Understanding schema locks in File Geodatabase registered to ArcGIS Server?

texts:
I have a folder containing many File GDBs. I registered this folder with my ArcGIS Server machine. I then published some mxd from ArcMap with data from some of those same fGDBs. The service is working fine and no data is copied to the server (registering worked fine).

I can see many new locks in my fGDBs, and I understand they come from the published active services.

As I can see from here, locks are of various types:

RD-Lock - Read Lock
SR-Lock - Schema Lock
WR-Lock - Write Lock
ED-Lock - Edit Lock

Mine are actually SR-Lock type. From ArcMap, I am able to start an edit session on some of the already published data, modify them, and I can see those changes in the service after saving, as expected.

However, I have troubles when it comes to (e.g.) adding a new field. Is this type of operation not allowed when there is a schema lock?

If so, would stopping the service be enough to release the schema lock, performing the changes, overwriting the data and activating it again?

The reason why I am asking is to understand if this is the right way to proceed, as last time I tried, I think I messed up something and I received an error during publishing (Staging the service succeeded, but uploading it did not).


texts: How do I make the Layer and Layers Order boxes of a usable size in QGIS.  They come up as very small boxes.  All other panels appear to work OK  I use a Mac running High Sierra


texts: Fixing displacement apparently originated from image extent

texts:
I am having problems with misalignment since. since I started to use GEE. I think I know what's wrong, but I don't know how to fix.

var amz18 = ee.FeatureCollection('users/trasegis/amazonia_pasture_2018')
var pasture = ee.Image();
for (var yr= 2001; yr <= 2018; yr = yr + 1) {
var pasture = pasture.addBands(ee.Image('users/vieiramesquita/PASTAGEM3/pasture_lapig_'+yr).remap([0,15],[0,1]).rename('pasture_'+yr));
}

Map.addLayer(pasture.select('pasture_2017'),{min:0,max:1},'raster2017') // example of a band without displacement
Map.addLayer(pasture.select('pasture_2018'),{min:0,max:1},'raster2018') // last band with displacement

Map.addLayer(table,{},'vector2018') //original vector for pasture_2018

https://code.earthengine.google.com/efd3db90c2abce545ad255efa45a2329

In the example above all the bands, except by the last one pasture_2018, are in the same crsTransform parameters. But the last one is different and clearly with a displacement (apparently 1 cell right and 1 to the top).

How can I fix this inside GEE? Is it possible?


texts: How to simulate google maps user experience

texts:
I’m not familiar so much with GIS applications and libraries. There is something that I would like to know if it’s possible or not.

We are using ESRI suite of products as our GIS infrastructure. We use products like ArcGIS desktop and ArcGIS server. However, the published maps are not that user friendly and not smooth. ESRI provides a lot of tools and utilities but not in the UX part. We want to reach a level when maps UX is something similar to the latest Google Maps.
Everyone can agree Google Maps has very smooth graphics, nice display to streets and amazing markers and info box animation.
I would like to ask if it’s possible to approach this goal or not by using some specific libraries or technologies.


texts: How to get a list of adjacent lots to a road in ArcMAP and PostGIS/QGIS?

texts:
In the following map, I have:

A line / road layer
A parcel / lot layer

The Task:

Find and count the list of lots adjacent to the road with their types (parcel_type is a field in the parcel layer)
The data is both stored as shp files and as a postgis tables

The Issues:

The line layer and the lot layer do not intersect i.e. the lot layer has the right of way  specified in it of the actual width of the road; whereas, the road layer is just a line layer.
a parcel might be adjacent to a road on more than one location. think of it as a large lot, with a lot within it. or lots with double frontages. (there is an example of such a lot in the picture below
This has to be automated on all roads of the study area (roughly 20,000 segments)

Image Sample


texts: How to break markers drawing which is in progress

texts:
The following leaflet codes draws a series of markers one by one. I don't want to wait until all points are drawn. How to break the drawing process before it ends and clear the screen?

<!DOCTYPE html>

  <title>Break Draw</title>
  <meta charset="utf-8" />
  <link rel="stylesheet" href="leaflet.css" />

  <link rel="stylesheet" href="Popcen.css" />

          <input type="button" id="btnSummary" value="Hit Sequence" onclick="myHit()"     class="btnStyle span2" />


texts: Keeping up to date with OSM changes in pgRouting

texts: I'm working on a routing application using OSM and pgRouting.  I'm only interested in a subset of data, mostly North America (performance was difficult for large amounts of data, so don't want to use the whole world).  I want to be able to keep up with updates to OSM data, but I haven't seen any place that provides a smaller changeset then for the whole world.  Is reimporting the data every time I want to update the only way?


texts: Error in opening Raster dataset in qgis "Cannot open GDAL Dataset"

texts: Error in opening Raster dataset in qgis (see image below). Is there any suggestion as to why this is happening? and how do i correct it?


texts: ArcObjects .NET - How to close/release FeatureClass, Workspace, Factory

texts:
I have a long-running process. I want to prevent resource leaks or rogue database connections.

At intervals during the process I want to do this:

get an ArcSDE (Oracle) workspace factory
open a workspace from the factory (at which point I get an open database connection)
get an existing feature class or table on the workspace,
query the feature class or table, loop over the cursor doing my business
then release/close everything such that:

The database connection and table lock from ArcSDE/Oracle's perspective (as revealed by something like "sdemon -o info -I users" or a query of the sde.table_locks table) is closed/released.
the process is resilient to ArcSDE/Oracle restarts (that is, I'm not leaving something hanging that won't work later after the nightly restart)
Any RCW, COM references, and memory are released.

Basically, due to the long-running nature of the process, I want to be really sure I haven't got resource leaks or rogue connections, and my process can survive ArcSDE/Oracle restarts.

I've seen discussions such as:

What are the rules for releasing ArcObjects from memory in .NET?
what every arcobjects programmer should know about singletons
how to release COM references
interacting with singleton objects

And this, from which I quote

  Each workspace factory maintains a pool of currently connected, active workspaces that are referenced by the application. When any of the previously listed Open* methods are called, the workspace factory verifies if a workspace has previously been opened with a matching set of properties. If so, a reference to the existing instance is returned.

All of which suggests to me that I should release (e.g. ComReleaser class or equivalent Marshal.ReleaseComObject() loop), probably in this order:

cursor
featureclass/table
workspace
workspace factory

Then there are discussions like this where people do all that, and maybe even sprinkle in System.GC.Collect() and their database connection still lives.

Oh gurus, what is the final straight dope on this?


texts: HTTP POST request: From cURL command to GeoEvent Input Service

texts:
I am obviously lacking some skills with the HTTP requests, yet I have to implement it into GeoEventServer (GES). The request with cURL is pretty simple, and does not require authentication:

curl -X POST https://<api-url>/search -d {} -H 'Content-Type: application/json;charset=utf-8'

It works fine in a command, and outputs a JSON string.

However, no output is created in GES:

What am I missing?


texts: Using multiple wildcard parameters with arcpy.ListFeatureClasses

texts:
I know there has to be a way to do this but I'm having trouble figuring out how to do so. I would like to create a list of feature classes to perform operations on however, I want to have the list contain feature classes which have 2 separate and distinct file name identifiers. The code would look something like:

FCS = arcpy.ListFeatureClasses("X_*" OR "*_Y")

I know OR is incorrect however, is there some way to do this without having to re write code to work with the second wildcard. Thanks for helping gout, i know there is a selection method in the Select By Attributes which allows for multiple wildcard use using the OR statement.


texts: How can I add button to the QGIS toolbar or create my own toolbar?

texts:
I want to create a button and add it to the QGIS toolbar or alternatively create my own toolbar and add the button here.

A plugin should be started when the user clicks the button.


texts: Active Directory AND integrated user store on ArcGIS Server 10.2

texts:
We have one License for ArcGIS Server 10.2.

We need to publish Services for internal use with Active Directory authentication and for external use (internet) with username/password authentication.

What would be the best way to achieve that?

Is it possible to create two sites with one AGS License on one server so that I can use different user-stores?


texts:
So after writing in a projectfile as in the cookbook in QGIS 2.14 with.

finfo = QFileInfo(filename)
project.write(finfo)

I want to open that project in the main window. So  I try:

project.read(finfo)

But the project doesn't open.
How can I open the .qgs file I've created in the window?


texts: Determine in Python if point is within a WGS84 polygon intersecting 180th meridian

texts:
I have a Python list of WGS84 points. Connecting these points together, we get two polygons on the WGS84 model of the Earth. I am interested in the smaller polygon which can possibly intersect 180th meridian and I would like to check if some given point is within it.

Please, note that I am already aware of the solutions that work well for the cases where a polygon does not intersect 180th meridian such as
Determine if point is within an irregular polygon using Python
or Determining if coordinate is within polygon?

I would like to get a library solution if possible so that I do not have to write additional tests for this function.


texts: Projection ETRS89 (Ellipsoid: GRS80) 6-Degree Zone 33 to UTM33N - WGS84

texts:
I'm having trouble to reproject a shape file in Qgis.

The projection is supposed to be "ETRS89 (Ellipsoid: GRS80) 6-Degree Zone 33" and I would like to project it to UTM33N - WGS84 (EPSG 32633).
When I set the Layer CRS to EPSG 25833 or EPSG 4258 and save it as UTM33N-WGS84 I don't get the result as intended...

This is the file I would like to reproject:
http://www.projektidee.org/images/extern/landschaftsschutz.zip

It should fit to this area in the end:
http://www.projektidee.org/images/extern/gebiet_UTM33.dbf.zip

Maybe one of you could give me a hint what I'm doing wrong?


texts: Tool or program to clean polygon dataset

texts:
I have a polygon dataset (shapefile, mapinfo tab file) which I would like to "clean" to be topologically correct. The polygons are assigned an ID which represents a specific class. For example ID 3 is forest and ID 2 is grass. In general the polygons does not touch or overlay but there are some exceptions. For example, a small polygon, forest, is on top of a larger polygon, grass. What I would like to do is to have a hole in the grass polygon where the forest polygon intersects.

What I have tried thus far was the v.clean tool in QGIS/grass plugin. But for some reason it also creates a grass polygon the same size as the forest polygon. So instead of two polygons there are three polygons. The images below illustrate what I have, what I want, what I would like and what the v.clean tool does.

Illustrates the tree polygon on top of the grass polygon

The forest polygon moved to illustrate that the grass polygon does not have a hole

The result from v.clean. Three polygons and a hole.

The result I would like

There are about 100 000 000 polygons (not all in one file) so doing it manually is not an option. I have arcmap and mapinfo but would like an opensource (QGIS/greass/saga/postgis/python script) to clean the polygons. Any suggestions?


texts: Why do Shapefiles and GeoJSON behave differently in GDAL Python?

texts:
I have a test.shp. I loop over the file, buffer the features, and write the buffers to a new file. I tested writing the buffers to a  Shapefile and to a GeoJSON. If I try to use the new buffer layer afterwards, for example to get a feature count, I get a number for the Shapefile 90 features, but not for the GeoJSON 0 features. What is the reason for that?

The code looks like this

import ogr
test = ogr.Open('test.shp', 0)
lyrTest = test.GetLayer()

For the Shapefile

shpdriver = ogr.GetDriverByName('ESRI Shapefile')
ds = shpdriver.CreateDataSource('Buffer.shp')
lyrBuffer = ds.CreateLayer('buffer', geom_type=ogr.wkbPolygon)
featureDefn = lyrBuffer.GetLayerDefn()

featureTest = lyrTest.GetNextFeature()
while featureTest:
  geomTest = featureTest.GetGeometryRef()
  geomBuffer = geomTest.Buffer(250)
  outFeature = ogr.Feature(featureDefn)
  outFeature.SetGeometry(geomBuffer)
  lyrBuffer.CreateFeature(outFeature)
  outFeature.Destroy()
  featureTest = lyrTest.GetNextFeature()

print lyrBuffer.GetFeatureCount()

For the GeoJSON

GeoJSONdriver = ogr.GetDriverByName('GeoJSON')
ds = GeoJSONdriver.CreateDataSource('Buffer.geojson')
lyrBuffer = ds.CreateLayer('Buffer.geojson', geom_type=ogr.wkbPolygon)
featureDefn = lyrBuffer.GetLayerDefn()

featureTest = lyrTest.GetNextFeature()
while featureTest:
  geomTest = featureTest.GetGeometryRef()
  geomBuffer = geomTest.Buffer(250)
  outFeature = ogr.Feature(featureDefn)
  outFeature.SetGeometry(geomBuffer)
  lyrBuffer.CreateFeature(outFeature)
  outFeature.Destroy()
  featureTest = lyrTest.GetNextFeature()

print lyrBuffer.GetFeatureCount()


texts: How to disconnect a feature from Geometric Network using Python ArcGIS

texts:
I am creating a Python add-in for ArcMap, which on click, will disconnect the feature from Geometric Network and then delete it from the feature class. But ArcPy doesn't have any method to disconnect the feature. So, I looked in ArcObjects and there they have method disconnect in INetworkFeature. I am a Python developer and don't know ArcObjects ecosystem. How can I do this?

Assumption of this add-in is that Geometric Network is added to ArcMap, the editor is started and you have selected the feature which you want to disconnect and delete from feature class.


texts: Cache FeatureCursor results so they can be requeried

texts:
What is the best method for caching features obtained by a featurecuror so they can be requeried without having to round trip to the datasource?

I am querying a subset of features using a queryfilter and later in my code want to subset the original subset of features (ie. I don't want to make a round trip as I already queried the features once). Is there an object or something that allows you to use queryfilter on a subset of client cached features? Or is just storing the features in a List of something and do searching via a loop?

In an arcuser article (fall 2009) http://www.esri.com/news/arcuser/1009/files/5tips.pdf
which has arcobjects optimizations there is a section on the following:

Code Listing 4: Don’t Be Too Selective (or Chatty)

which offers:

Build a where clause containing all items you want to fetch.

Create a query filter and set the where clause property to the
where clause you just created.
Create a Geodatabase RecordSet and set the source table
property.
Create a new InMemoryWorkspace and save the RecordSet
as a table within it.

but not much more than that. I can't track down the code samples for the article and the documentation on IRecordSet (assumption that is what this article is referring to) is very light and not much to go on.


texts: Solve for alternate routes without barriers or additional stops in ArcGIS Network Analyst?

texts: I created a network dataset using ArcGIS Network Analyst. I have one start and one destination point. I want to find the three shortest paths from start to my destination (find three shortest paths from point a to point b). Is it possible to run an analysis for three different routes?


texts: Set HTML ID or class for an icon in OpenLayers

texts: I wanted to know if there is any way to set an HTML ID for an icon in OpenLayers. For example, If I set an icon with HTML ID "icon1", I want to be able to access it using document.getElementById("icon1");


texts: What is a good Javascript web API (preferably open) for building web maps that require little/no maintenance?

texts:
I have recently volunteered to make a web map for a local non-profit. The map will be pretty basic. It needs to display Lane county, Oregon split into multiple districts, with a "central office" location for each district. I plan on building this in my spare time over the summer, and I want to be able to hand this off to them with very little risk of it needing any sort of future maintenance.

I have built a number of web maps using the ArcGIS Javascript API with ArcGIS Server, but I am planning on learning a different API since this organization doesn't have any ESRI software (or any GIS for that matter). I have been considering OpenLayers/OpenGeo. One thing to keep in mind is that this organization doesn't/won't have a server set up to host the layers of this map so I'm not sure how I would store and access the data. JSON maybe?

I'm currently unfamiliar with most non-ArcGIS solutions, so any advice would be appreciated.


texts: Get pixel(Row Column) from Lat Long for Moon

texts:
How to obtain pixel(Row Column) Value from Latitude and Longitude given, for Moon?

Actually, I am Having a GeoTiff file which contains craters, and I need to crop the specific Crater whose Center is given to me in the form of Latitude and longitude.

I know that I have to use gdal to get this done, but I am having wrong results from my code.

Code is in Python
Code:

from osgeo import osr, gdal
import sys

# get the existing coordinate system
#SPHEROID["MOON_localRadius",1737400,0]
# create the new coordinate system
ds = gdal.Open("C:\Users\ZARL\Desktop\PRL_Project\small.tif")
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())
wgs84_wkt = """
        GEOGCS["GCS_MOON",
    DATUM["D_MOON",
        SPHEROID["MOON_localRadius",1737400,0]],
    PRIMEM["Reference_Meridian",0],
    UNIT["degree",0.0174532925199433]],
PROJECTION["Equirectangular"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",22.46049118041992],
PARAMETER["standard_parallel_1",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs)

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5]

#get the coordinates in lat long
#latlong = transform.TransformPoint(minx,miny)

transform1 = ds.GetGeoTransform()

xOrigin = transform1[0]
yOrigin = transform1[3]
pixelWidth = transform1[1]
pixelHeight = -transform1[5]

latlong = transform.TransformPoint(xOrigin,yOrigin)
#data = band.ReadAsArray(0, 0, width, height)

points_list = [ (22.4655,-71.7136)] #list of Lat long

for point in points_list:
    col = int((point[0] - latlong[1]) / pixelWidth)
    row = int((latlong[0]-point[1] ) / pixelHeight)

    print row,col

I am getting wrong values.

The Tiff file is:https://www.dropbox.com/s/twov5jpe40hcf4q/small.tif?dl=0

when I run gdalinfo command on the picture I get the following:

    Driver: GTiff/GeoTIFF
Files: small.tif
       small.tif.aux.xml
Size is 40, 33
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["GCS_MOON",
        DATUM["D_MOON",
            SPHEROID["MOON_localRadius",1737400,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",22.46049118041992],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (135.724395353784870,-2174581.951404873300000)
Pixel Size = (0.812720930262190,-0.812720930262190)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (     135.724,-2174581.951) ( 22d27'53.88"E, 71d42'47.22"S)
Lower Left  (     135.724,-2174608.771) ( 22d27'53.88"E, 71d42'50.40"S)
Upper Right (     168.233,-2174581.951) ( 22d27'57.74"E, 71d42'47.22"S)
Lower Right (     168.233,-2174608.771) ( 22d27'57.74"E, 71d42'50.40"S)
Center      (     151.979,-2174595.361) ( 22d27'55.81"E, 71d42'48.81"S)
Band 1 Block=40x1 Type=Byte, ColorInterp=Gray
  Min=1.000 Max=254.000
  Minimum=1.000, Maximum=254.000, Mean=122.977, StdDev=48.776
  Metadata:
    STATISTICS_MAXIMUM=254
    STATISTICS_MEAN=122.97651515152
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=48.776222721807

Any Other Approach is also fine.


texts: How to integrate multiple tile outputs generated from gdal2tiles into single google map?

texts:
Todo

Generated tiles from image1.tif - done
Generated tiles from image2.tif - done
Need to overlay these 2 sets of tiles on top of google maps - this is the problem

gdal2tiles generated ready-made HTML that plotted those generated tiles but since it came from 2 different image files, it also generated 2 ready-made HTMLs.

How do I integrate those multiple tile outputs generated from gdal2tiles into a single google map?


texts: Collapse points onto route, get distance of each point from start of route

texts:
I have a polyline SHP (with only one feature) that I have converted to a route using arcpy.CreateRoutes_lr(). I have a separate points SHP that is spatially scattered about the line. I want to "collapse" the points onto the route, and then for each FID in the points file, find the distance of the collapsed point along my route, relative to the start of the route.

I tried Locate Features Along Routes just in ArcMap from the toolbox, as suggested by "gisnside" here, but the distance returned is the distance between the point and the line (exactly the same output as Near Analysis gives), not the distance of the collapsed point along the line.

Other details:

Both polyline and points SHPs are in the same UTM projection.
The "route_id" layer is correct in the route file.
All of the points have a unique integer id in a "pointid" field
Must be accomplished in ArcPy; process to run many times.

# file location setup
myline = r'mypath\line_projected.shp'
mypoints = r'mypath\points_prpjected.shp'

arcpy.AddGeometryAttributes_management(myline, "LENGTH_GEODESIC", Length_Unit = "METERS")
outroute = mypath + '\\' + 'myroute.shp'
arcpy.CreateRoutes_lr(myline, "route_id", outroute, "LENGTH", coordinate_priority = "LOWER_LEFT")
# ... what tool to use ...


texts:
I want to deploy multiple GeoServer instances using AWS Elastic BeanStalk. In order to do this, I figured I have to keep my data dir somewhere else, like S3. I suppose the naive way to do this is to set up a job to run "s3 sync" every 30 seconds or something like that.

Has anyone ever solved this problem?


texts:
I cannot solve the mystery of no outputs in any case and in any experiment I have done with all kinds of parameters in the latest QGIS 3.2.3 Bonn and GRASS 7. r.texture does not provide any output. Not in temp folders, not in user defined folders. I have not used this function before to know the details however.

Anyone else having the same issue? Any workarounds?

The last time I run it, I got the following LOG output:

  Processing algorithm… Algorithm 'r.texture' starting… Input
  parameters: { '-a' : False, '-s' : False,
  'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' :
  None, 'distance' : 1, 'input' :
  'C:/Users/.../AppData/Local/Temp/processing_b31ae406d97346319858b8acfa3fe9ea/d40ba480206c44e1855a50ce4442b154/output.tif',
  'method' : [0], 'output' : 'C:\Users\...\textures', 'size' : 3 }

  g.proj -c proj4="+proj=utm +zone=...+datum=WGS84 +units=m +no_defs"
  r.external
  input="C:\Users...\AppData\Local\Temp\processing_b31ae406d97346319858b8acfa3fe9ea\d40ba480206c44e1855a50ce4442b154\output.tif"
  band=1 output="rast_5bc7529c4be7121" --overwrite -o g.region
  n=...s=...e=...w=...
  res=0.08917000000001485 r.texture input=rast_5bc7529c4be7121
  method="asm" size=3 distance=1
  output=output0721b2f6bc07456a8304b894c5bcf6e4 --overwrite for r in
  $(g.list type=rast pattern='output0721b2f6bc07456a8304b894c5bcf6e4*');
  do r.out.gdal -m -t input=${r}
  output=C:\Users...\textures/${r}.tif
  --overwrite -c createopt="TFW=YES,COMPRESS=LZW" done Starting GRASS GIS...

  WARNING: Concurrent mapset locking is not supported on Windows

  Executing

  ...

  C:\OSGEO4~1\bin>chcp 1252 1>NUL

  C:\OSGEO4~1\bin>g.proj -c proj4="+proj=utm +zone=...+datum=WGS84
  +units=m +no_defs"

  Default region was updated to the new projection, but if you have
  multiple mapsets g.region -d should be run in each to update the
  region from the default

  Projection information updated

  C:\OSGEO4~1\bin>r.external
  input="C:\Users...\AppData\Local\Temp\processing_b31ae406d97346319858b8acfa3fe9ea\d40ba480206c44e1855a50ce4442b154\output.tif"
  band=1 output="rast_5bc7529c4be7121" --overwrite -o

  Over-riding projection check

  Reading band 1 of 1...

  r.external complete. Link to raster map
  created.

  C:\OSGEO4~1\bin>g.region n=...s=...
  e=...w=... res=0.08917000000001485

  C:\OSGEO4~1\bin>r.texture input=rast_5bc7529c4be7121 method="asm"
  size=3 distance=1 output=output0721b2f6bc07456a8304b894c5bcf6e4
  --overwrite

  Reading raster map...

  Calculating Angular Second Moment...

  2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

  r was unexpected at this time.

  C:\OSGEO4~1\bin>for r in $(g.list type=rast
  pattern='output0721b2f6bc07456a8304b894c5bcf6e4*'); do

  Execution of

  finished.

  Cleaning up temporary files...

  Press any key to continue . . .

  Starting GRASS GIS...

  WARNING: Concurrent mapset locking is not supported on Windows

  Executing

  ...

  C:\OSGEO4~1\bin>chcp 1252 1>NUL

  r was unexpected at this time.

  C:\OSGEO4~1\bin>for r in $(g.list type=rast
  pattern='output0721b2f6bc07456a8304b894c5bcf6e4*'); do

  Execution of

  finished.

  Cleaning up temporary files...

  Press any key to continue . . .

  Execution completed in 7.08 seconds Results: {'output':
  'C:\Users\...\textures'}

  Loading resulting layers Algorithm 'r.texture' finished


texts: Raster vs. Raster Scatter Plot in ArcGIS

texts:
I created two raster datasets from the same world view 2 image. One dataset is NDVI, the other is a normalized surface water index (NDSWI). I would like to create a scatter plot for a pixel by pixel comparison.

I've seen a few options for how to do this but they involve running a combination function in raster calculator which only allows for integers. I need the float data. I also have access to ENVI if there is a way to do it in that software.

Running ArcGIS 10.2 and ENVI 4.4


texts: Looking for (free) WMS of surface waters of Europe

texts:
Do you know any WMS with hydrographic map of Europe?
Especially I am looking for map with surface waters of Europe


texts:
I'm trying to completely remove and reinstall QGIS on my system (Mac os 10.11.6). I follow the steps in the Kyngchaos and delete the existing application. But on reinstalling framework, components and program I still see old plugins from prior installation. I can't delete these which i what i want to do. As screenshot, they are there but i can't uninstall or reinstall as both boxes are greyed out. I've installed the latest Kyngchaos framwork and build 2.18.15. Please can someone advise where these are stored on my system so i can delte them.

EDIT by looking in the package contents and going to folder Plugins/qgis I can remove all but these which aren't in that folder. If i delete everything in the folder it causes QGIS to break and I have to reinstall;


texts: Automatic method for placing values of determined column from a SpatialLinesDataFrame near contour lines

texts:
I have following vector layer that represents countour lines. With QGIS 3, labels placed near lines (Layer Properties -> Labels) were taken from TEMP_ENE field in attributes table.

I would like a similar behavior in R. In this case, I used rgdal library to read vector layer as follows:

library(rgdal)
isoterma.ene<- readOGR('isotermas_1grado_ene_4326.shp')
isoterma.ene
class       : SpatialLinesDataFrame
features    : 24
extent      : -71.84042, -68.01765, -41.09749, -36.15523  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 2
names       : ID, TEMP_ENE
min values  :  0,       12
max values  :  8,       23

I found out this question where is used locator() function.

coords <- locator()
coords
$`x`
[1] -68.97818
$y
[1] -38.4115
text(x=-68.97818, y=-38.4115,label="23")

I tried it out in only one case (value 23) and it worked.

Is there one method where automatically values in isoterma.ene$TEMP_ENE of SpatialLinesDataFrame can be placed in these lines?

Editing Note:

Based in mdsummer answer, following approach works manually:

library(rgdal)
library(raster)
isoterma.ene<- readOGR('isotermas_1grado_ene_4326.shp')
plot(isoterma.ene)

By using click from raster library:

c1<- click(isoterma.ene, id = TRUE, xy = TRUE)
c2<- click(isoterma.ene, id = TRUE, xy = TRUE)
c3<- click(isoterma.ene, id = TRUE, xy = TRUE)
values<- c(c1$TEMP_ENE, rev(c2$TEMP_ENE), rev(c3$TEMP_ENE))
values
[1] 23 22 21 20 19 18 17 16 15 14 13 12 14 13

It looks as follows:

By using locator():

coor1<- locator()
coor2<- locator()
coor3<- locator()

x<- c(coor1$x, coor2$x, coor3$x)
y<- c(coor1$y, coor2$y, coor3$y)
x
 [1] -69.04728 -69.47572 -69.62774 -69.71067 -69.86269 -70.00090 -70.15292
 [8] -70.24967 -70.38787 -70.66428 -70.89923 -71.12036 -70.12528 -71.17564
y
 [1] -38.50867 -39.16729 -39.53439 -39.73954 -39.86910 -39.96628 -40.13903
 [8] -40.33338 -40.50613 -40.61410 -40.75446 -41.09997 -36.26288 -37.01868
plot(isoterma.ene)
text(x=x, y=y,label=values, cex=0.6)

Result is similar as expected; but in manual way:


texts: Seeking shapefile of administrative divisions of Cote d'Ivoire (Ivory Coast) post-2011

texts: I've been trying to find a shapefile containing CIV's administrative divisions (preferably down to the level of sous-prefecture) reflecting the national reorganization of 2011 but have had no luck. Everything I've found thus far contains pre-2011 boundaries; does anyone know where I might find a shapefile representative of 2011 boundaries?


texts: Use Parse Path and Get Field Value for Output FC Name? Model Builder

texts:
Is it possible to combine my parse path "value" and my get field value "date" to name my output feature class?

Essentially, I want the name to appear as such: "_40588_04102016"
So my output feature name is: %Value%_%Date%

When the tool is run, the output name appears as such: "_40588_1"

Can this be done?
If so, what am I doing wrong?


texts: Extract the "mid/center-point" of the touching borders/lines of adjacent polygons

texts:
I have a layer with several polygons and want to extract the "mid/center-point" of the touching borders/lines of adjacent polygons. Despite trying several methods, I do not know how to achieve what I had wanted.

What I tried recently:

I used Polygons to lines to convert the polygons to lines
I then used Explode lines to split the lines at vertices

But Explode lines just splits at all of its vertices. I am thinking that the steps might be to:

Convert the polygons to lines
Split the lines at vertices only if there is one or more line touching the vertex (like splitting the lines only where the black dots are as shown in the screenshot below):

Find the mid-point (point coordinates) of each of the resulting split lines (which would be the red dots below):

What are the steps and algorithms needed to achieve this (extract/get the red dots as shown in the image)?

Update: I tried QChainage plugin but not all common vertices were extracted (like these green circles):


texts: Use Geokettle to Transform file geodatabase data

texts:
I would like to use GeoKettle to transform my file geodatabase data into something else (most likely postgis).  Admittedly I am very new to ETL but I am absolutely stumped when I try to define an input as file geodatabase.  "ogrinfo --formats" list file geodatabase as an option.

The reason I want to use the ETL is my long term solution will have to do with some schema changes (field names and probably types ) and with that mapping attributes to the new layer.

All of that is fine and good except I can't get out of the gate with using a file geodatabase as my input.

Ubuntu 14.04

Geokettle 2.5


texts: Error updating geopackage layer attribute in Python using sqlite3

texts:
I'm attempting to update an attribute field in a geopackage layer using the Python sqlite3 module, but any attempt at update is giving an error "no such function: ST_IsEmpty".
Is there anything I can do to make the updates work?
Running python-qgis in an OSGeo4W shell (python 3.7):
import sqlite3
gpkg = 'addressbaseplus.gpkg'
conn = sqlite3.connect(gpkg)
cur = conn.cursor()
SQL="update os_addressbaseplus_joined set address = 'xxxx' where fid = 22"
cur.execute(SQL)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
sqlite3.OperationalError: no such function: ST_IsEmpty

Google tells me this error is something to do with spatial indexing, but I'm not going near the geometry column???


texts: Using option "resize" when making temporary layer into permanent shapefile in QGIS

texts:
I have a temporary layer I want to make permanent into a .shp file. In the settings, there is an option called resize. I can choose between yes and no.

What is the usecase of this option and when should I use "Yes"?


texts: What level of error expected when using NAD StatePlane CRS from adjacent state?

texts:
I am adding Georgia data to my GIS database and, unlike the California data I've added, have decided to use my existing EPSG:3361 (South Carolina StatePlane) coordinate reference system for calculations.

I've done this with reasonable success in North Carolina (South Carolina's northerly neighbor) but never did the research to determine the error suffered by using a linear reference system from a non-intended, albeit adjacent area.

I'm guessing the error shouldn't be more than 50 feet, am I correct in assuming this?

More technical information:

here's the linear reference system I'm using in Georgia: http://spatialreference.org/ref/epsg/3361/
here's the linear reference system I should be using in Georgia but am not: http://www.spatialreference.org/ref/epsg/2239/


texts: MapInfo Thematic Mapping - All other values?

texts: Relatively quick question here, how do I go about adding a style for all other values using individual value thematic mapping in MapInfo?


texts: Does arcpy.da.SearchCursor store the results in an array?

texts:
When performing a query like this:

var = 2

with arcpy.da.SearchCursor(in_table="table", field_names=['c1','c2'], where_clause = '[c1]' + '=' '\'' + str(var) +'\'') as cursor:

I've just simply been doing a for loop below this to do a bit of value checking.

var = 2

with arcpy.da.SearchCursor(in_table="table", field_names=['c1','c2'], where_clause = '[c1]' + '=' '\'' + str(var) +'\'') as cursor:
     for x in cursor:
         if x[1] != 1:
              print whatever

Is it possible to eliminate the for loop by doing something like:

  if cursor[0][1] !=1:
     print whatever

I may not be understanding this fully, but is this query not storing the columns that I have specified as an array into the variable cursor? And should I not be able to reference the second column of the query utilizing cursor[0][1]?

I've tested this out and it doesn't work, so I'm guessing that it's not doing what I'm expecting to do so, I just do not know how to access the result without doing a for-loop.

Sorry if this question is confusing.


texts:
Using ArcPy, how do I sum a range of fields?
While I know this is possible using field calculator, I am looking for something repeatable for long list of fields that changes. I already have code that will determine which fields should be included and gets placed in a list with one additional field that will house the total. This list is the variable workfields The idea of what I'm looking for looks something like this:
with arcpy.da.UpdateCursor(fc, workfields) as cursor:

    for row in cursor:
        row[-1] = sum(row[0]:row[-2])
        cursor.updateRow(row)

del row, cursor

Basically row[-1] = sum(row[0],row[-2]) needs to work like:row[-1] = row[0] + row[1] + row[2] + .. row[-3]+row[-2]. I have gotten come code to work that loops through the entire list updating a variable using y += row[i] that can then be passed through as to what row[-1] should be. I think that going this route is a bit clunky and just adds too much time in running the script when I have hundreds of fields I need added up for thousands of records.


texts: Extracting Address numbers in a Block Reference from autocad

texts:
I have an autoCad file of a subdivision and I want to extract the addresses from it. However, when I export the layer to a SHP, it loses any information it had about the actual addresses. The addresses are stored as a Block Reference Layer.

Usually, I can use Explode or Burst to break up the layer and then export. However, those two options are making AutoCad crash instead.

Does anyone have any recommendations on other ways to obtain the data? I'd rather not enter it all manually into ArcMap.

Edit:
Solutions Tried

1) I tried using the Import Cad Annotation tool and got an Error 000278 and 000016. This looked to be, maybe, due to an extent error which I've seen before when using the generated GDB's we get for these jobs. I made a new GDB just for the cad file to try to work around the extent issue and re-ran the tool. ...Only to get the same two errors. I am allowing the tool to create its own feature class but the error still occurs.

2) I added the annotation via Add Data and then right clicked and using Import Cad Annotation. This produced the same two errors.

3) In AutoCad 2012, I used WBLOCK to create a new cad file with ONLY the layer I wanted (Address). This ended up changing the data from an annotation layer to a CAD point file. Copy CAD Feature layer. This finally produced the closest result to what I wanted, a feature class with an attribute that had the address numbers in it.


texts: ArcGIS Pro Field Calculator completing without error, but no change to Feature?

texts:
I'm attempting a research project in which I am seeking to categorize a municipality's Property Use Codes as Generalized Land Uses.  For example, "AGRII - Developed Farmland" and "INDGR - Grain Elevator" are both "Agricultural" land uses.  I have created a new field "Land Use Category" and used the following code:

PLY_AssessmentParcels_LandUseCategory = MyFunction(!Assessment_Parcels_CSV_csv_Property_Use_Code!)
_

def MyFunction(s):
   if ("AGRII" or "INDGR") in s:
      return "Agricultural"
   elif ("VRES2" or "VRES1" or "VINDU" or "VCOMM" or "VAPRK" or "VAGRI" or "CNVAC") in s:
      return "Vacant"
   elif ("STATR" or "SKYWK" or "RAILR" or "ARPRT" or "CMPST" or "CMPSP") in s:
      return "Transportation"
   elif ("RESTR" or "RESSU" or "RESSS" or "RESSD" or "RESRM" or "RESRH" or "RESOT" or "RESMU" or "RESMH" or "RESMC" or "RESMB" or "RESGC" or "RESAP" or "RESDU" or "RESAM" or "CNAPT" or "CNCMP" or "CNDRH" or "CNRES") in s:
      return "Residential"
   elif ("CMCMU" or "CMFBK" or "CMHBH" or "CMHHO" or "CMHMO" or "CMMRH" or "CMOFF" or "CMOGV" or "CMOMC" or "CMRCS" or "CMRCV" or "CMRNS" or "CMRRE" or "CMRRS" or "CMRSM" or "SMRST" or "CMRWC" or "CMVCD" or "CMVSR") in s:
      return "Commercial"
   elif ("CNCOM" or "REFRL" or "CMSTP") in s:
      return "Other"
   elif ("PIRST" or "PIRRT" or "PIRRI" or "PIRPO" or "PIRPK" or "PIRMU" or "PIRGC" or "PIRCC" or "PIRBA" or "PIRAR" or "PIEBM" or "PIECS" or "PIETC") in s:
      return "Recreational"
   elif ("CNIND" or "HYDSS" or "INMHM" or "INMLM" or "INMMI" or "INMMU" or "INWSC" or "INWWH" or "MTSSS" or "TRAPI") in s:
      return "Industrial"
   elif ("PIEMA" or "PIICH" or "PIIDC" or "PIIGC" or "PIIHO" or "PIIMB" or "PIINH" or "PIIPF" or "PIIRE" or "PIISC" or "PIIUC" or "PIMCE" or "PIMFH" or "PIMLI") in s:
      return "Institutional"
   else:
      return

When I run the above code, I often received an error stating "Cannot acquire lock".  After restarting ArcGIS Pro, I can usually get one good run out in which the process completes successfully in about 40 seconds.  The problem is that nothing has changed in my attribute table.  Please note, I'm dealing with about 212,000 records.  The file is stored in a .gdb and I have re-created the feature class in the .gdb with the same issues.  System specs are very reasonable (i7 7800k, 32gb ram, GTX 1080i, etc.)

Any ideas?


texts: Determining closest representation of polygon as collection of points+radii

texts:
Have a bit of a tough question, which I realize there may not be an easy solution to. Figured what better place to ask.

I'm attempting to transform a polygon into a series of arbitrary points with arbitrary radii (circles!) that best represent the area covered by the polygon.

Minor over- and underflow is acceptable, as is circle-overlap, in hopes of achieving an efficient solution (i.e. fewest possible points, no gaps).

Ideally, a given poly would be represented by a few large circles, and several smaller circle on the perimeters.

Essentially, the problem is that I have a dynamic number of polys that get hit on geo-spatial queries given a specific gps coordinate, however, we are required to move to a system wherein I will not be able to utilize a point-within-poly query, but will have to rely on point-within-distance queries.

Hopefully someone has at least attempted something similar, and, if not, hopefully someone is willing to throw some ideas around!

Open to most languages, but this needs to be done programmatically!

Update per clarification in comments:

My point is that I feed indexes to the system and I get alerted when they are hit. I don't perform any queries on this system myself (black box), so I wouldn't have control enough to negate a query. That's the reason I need to transform the poly to a representation of points.


texts:
I have a number of LIDAR files (the extension is .adf). When put together, they show the same area (each shows one part of it), so, before I do anything else with them, I would like to combine them into one.

In other words, the equivalent of the Append operation on vector data. What's the best way to go about this in ArcGIS 10?


texts: How to align TIF and GPX layers in QGIS?

texts:
My initial layer is an OS tif file, with its associated tfw.  I can zoom into the area I want, no problem.  The CRS is OS 1936 (and that's the project CRS as well).

I download a 'track' layer from an Etrex as a GPX file: I make sure the CRS is OS 1936.  If I select this and zoom to the layer extent I can see it.

But I cannot see both layers at once (that is, overlaid) - although I achieved this with another similar dataset yesterday.

I'm obviously doing something wrong which I suspect is really very basic - can anyone tell me what?


texts: Remove layer switcher controller in leaflet

texts:
I have created a map in leaflet and added an overlay as:

map.layerscontrol.addOverlay(tile_layer,"{{ resource.title }}");
map.addLayer(tile_layer);

The map shows great and the layer as well.
Now I want to remove the 'layer switcher' controller.
How can I do that in Leaflet?


texts:
There is an analysis about the potential of solar areas in Germany. When I finished the project I realised, that there is one polygon with a huge value for its area. By looking closer it gets clear that this one geometry includes nearly 100 "splinters" all over the study area.

What I want to do is to cut off and remove all superfluous parts/splinters. It would be best if its possible to sort out by area of the splinters.


texts: Change in altitude value after resizing .tif file

texts:
I have a TIFF file of 7378x6923 (widthXheight) resolution. I resized it to 1200x1200 pixels using Image module in Python.

img = Image.open('path to tif file')
img = img.resize((1200,1200))

But after resizing, the pixel values show unusual values. Before resizing, a certain pixel's value was 125.28. Now, after resizing, that pixel's value has turned to 1123715645. I am reading the pixel value with:

with rasterio.open('path to tif file') as dataset:
    band = dataset.read(1)
    print(type(band))
    value_of_interest = band[row_of_interest, column_of_interest]
    print(value_of_interest)

Why does this happen?

How can I preserve the height values in same pixel position?


texts: How to add a polygon from PostGIS ST_AsGeoJSON into Leaflet?

texts:
I try to use the following function to add polygons to my Leaflet map.

function polygon(geojson)
{
    alert(geojson); // alert for testing
    layer = L.geoJson(geojson);
    layer.addTo(map);
}

When I call the above function with the GeoJSON string I fetched from a PostGIS table using ST_AsGeoJSON(), the following string is alterted.

{"type":"Polygon","coordinates":[[[10,370],[20,380],[30,410],[-10,395],[10,370]]]}

But no polygon shows up on the map. As far as I can say, there is no error message. And I can definitely say that a polygon with the points 10 370, 20 380, 30 410, 10 395 shows up on the map, if I insert it using L.polygon.

What's wrong?


texts: Viewing 2D subsurface seismic data (raster file) in QGIS using 2threejs 3d plugin

texts: I am an exploration geologist and trying to figure out if there is a way to view SUBSURFACE data using the QGIS 2threejs plugin? I have a seismic raster which I would like to view in 3D, with vectorized surface geology draped over. This is the kind of thing you would complete in leapfrog, but I would love to hear if anyone knows how to do this in QGIS!


texts: GeoServer, getting layer time dimension values

texts:
I want to consume a GeoServer layer that has the time dimension enabled from an external WMS service that I have no control over. Time dimension is updated frequently, several times in an hour. Using Leaflet to visualize this layer but I'm having trouble programmatically figuring out what are the values for time dimension at any point in time.
In the GetCapabilites request I see for example that
<Dimension name="time" units="ISO8601" multipleValues="0" current="1">2021-01-12T02:50:00Z/2021-01-12T09:40:00Z/PT5M</Dimension>
I'm wondering if there's a way to get GeoServer return only the possible time dimension values to the client. I could of course figure out some heuristics or parse the whole GetCapabilites request but this seems like killing a fly with nuclear bomb so looking for more elegant solution. Ideally an API endpoint that just returns an array of time dimension values for a layer.


texts: Combining Get Feature requests within OpenLayers3 Overlay

texts:
I have a map which currently shows one feature of a layer in an overlay popup when you click on it, as per this JSFiddle.

As you can see, it currently displays the address within the overlay. I am looking to show the 'address' and 'name' features within the same overlay, presumably by adding in feature requests.

Also, if possible this will go into a table within the overlay? However I'm not sure where I enter the code for a  to go into the overlay popup.


texts: How to write attribute-only tables in Geopackages with GeoPandas

texts:
I searched around and found similar questions (Writing Non Spatial Table into Geopackage - Python), but that discussion is on a slightly different topic.
Is there a way to write attribute-only tables into a Geopackage using GeoPandas?
You need to specify a geometry field to create a GeoDataFrame. This is the error you get if there isn't one:
Import geopandas as gpd

sql = "select * from attronly"
con = ###your connection###
gdf = gpd.read_postgis(sql, con)

ValueError: Query missing geometry column 'geom'

But you can create an empty geom field and then drop it afterwards. If you do that, however, there's a different error as the function still expects there to be geometries to parse:
Import geopandas as gpd

sql = "select null as geom,* from attronly"
con = ###your connection###
gdf = gpd.read_postgis(sql, con)
gdf.drop(columns=['geom'],inplace=True)
gdf.to_file("file.gpkg", layer='attronly', driver="GPKG")

AttributeError: No geometry data set yet (expected in column 'geom'.)

I suspect the solution is either using Pandas to write into a Geopackage (though I can't find any documentation about that), or going back and using Fiona directly. I can also write in that nulled geometry field and then remove it afterwards from the Geopackage, but I'm hoping for some expert advice on doing this without bandaid solutions.


texts: Getting coordinate system name from SpatialReference using GDAL/Python

texts:
In Python, using GDAL, I've extracted a raster's projection as a WKT string as follows:

wkt = dataset.GetProjection()
# wkt is 'PROJCS["GDA_1994_Transverse_Mercator",GEOGCS["GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.2572221010002],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]'

Using the WKT string, I can create a SpatialReference instance as follows:

src = osr.SpatialReference()
src.ImportFromWkt(wkt)

This is easy-peasy. I can extract several parameters of the projection from src such as the UTM zone, etc., easily enough. But I can't figure out how to extract something like the name of the projection, i.e. "GDA_1994_Transverse_Mercator". This has surely got to be possible, but the Python API documentation may as well be nonexistent for all the use it is.

How do I extract the names of the projection and the geographic coordinate system?


texts:
I would like to get a symbology using several columns. I have read that it was not possible to set up several fields within the CLASSITEM option in the mapfile. I have tried to do this in a different way but the layer is not displaying at all...

I'm using Mapserver 6.2.1. Here is a simplified extract of my mapfile's code :

CLASS
  NAME "Symbol 1"
  EXPRESSION (("[field1]"="01") && ("[field2]"="01"))
  STYLE
    SYMBOL circle
    COLOR 255 0 0
    SIZE 10
  END
END

CLASS
  NAME "Symbol 2"
  EXPRESSION (("[field1]"="02") && ("[field2]"="02"))
  STYLE
    SYMBOL circle
    COLOR 0 255 0
    SIZE 10
  END
END


texts: OD Matrix algorithm in QNEAT3 plugin (QGIS)

texts:
I have two questions regarding OD Matrix algorithms in QNEAT3 plugin.

If I set the cost strategy as the fastest path, would the output be the travel time? If yes, is the output in seconds?

The output of OD Matrix with the shortest path cost strategy is in meters?


texts: Standalone PyQGIS script is not working on Python virtual environment

texts:
Error while executing the script in Python virtual environment
(PVirtual) xor@HP-240:~/code/qgis$ python3 aa.py
Traceback (most recent call last):
  File "aa.py", line 2, in <module>
    from qgis.core import *
ModuleNotFoundError: No module named 'qgis'

working on terminal: if not using the Python virtual environment
xor@HP-240:~/code/qgis$ python3 aa.py
Warning: QT_DEVICE_PIXEL_RATIO is deprecated. Instead use:
   QT_AUTO_SCREEN_SCALE_FACTOR to enable platform plugin controlled per-screen factors.
   QT_SCREEN_SCALE_FACTORS to set per-screen factors.
   QT_SCALE_FACTOR to set the application global scale factor.
{'Pointmap': 0}
{'Pointmap': 5}
0 5
Point:  Pointmap Asset1 77.62277551171958 12.91716831378434

Since I'm using the MySQL and the Rest API also in the Python script so I need to use the Python virtual environment. but the QGIS code is not working in the Python virtual environment.