texts: How do I create a projection file using GDAL

texts:
I am trying to obtain a buffer from a polyline file using GDAL with the following function

def createBuffer(inputfn, outputBufferfn, bufferDist):
    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()
    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()
    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)
        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)

It works fine and creates the .bdf .shp .shx files

Well. But I would like to include in the same directory ( created by the function) the projectioin file .prj

I have got the following function.

def construye_proyeccion(nombre_arhivo):
    spatialRef = osr.SpatialReference()
    spatialRef.ImportFromEPSG(4326)
    spatialRef.MorphToESRI()
    file = open('%s.prj' % filename, 'w')
    file.write(spatialRef.ExportToWkt())
    file.close()

What can I do to create the buffer and his .prj file at the same time an in the same directory?


texts:
Using as an example Landsat 8 bands:

Does pansharpening band enhace only visible range and therefore can I use it only for analysis connected with visible spectrum? Or I can pansharpen every possible combination of raster composition so I can get an image with NIR/SWIR/Coastal/Cirrus included? If the latter is correct, will the algorythm perform some not-so-accurate 30->15m interpolation or it's 'pure' 15m? I see in the specs of L8 that Pan band ranges from 0.50 to 0.68 µm what covers G and R band what leads me to the first hypothesis, but I have seen some other compositions that were pansharpened.


texts:
I have a large point dataset of about 5.3 million points. I would like to delete the "lonely" points using QGIS, for example all points with no neighbor within a 5 km radius. The extent of my point dataset limits some of the tools I guess...

I tried the Distance matrix en v.neighbours tools, but QGIS can't handle it and goes irresponsive... Is there another way to work around this problem for large datasets?


texts: Showing size of pie chart in legend of ArcMap 10.7.1

texts:
I'm using a pie chart map in which the size of the pie chart convey an attribute value (scaled size).
Is it possible to show this classification in the legend of ArcMap 10.7.1 ?
There is a similar answer for QGIS but I can't 'adapt it to ArcMap 10.7.1 ( Showing size of pie chart in legend of QGIS? )
This is an example of a map with pie charts that I need to conclude them to the legend of the map.


texts: Extract field's names in R without reading the file

texts:
I was wondering if it's possible to use the R package sf for extracting the names of the fields in a file without reading it. For example, I would like to replicate the following behaviour:
system(paste0("ogrinfo -so ", system.file("shape/nc.shp", package = "sf"), " nc"))
#> INFO: Open of `/home/andrea/R/x86_64-pc-linux-gnu-library/3.6/sf/shape/nc.shp'
#>       using driver `ESRI Shapefile' successful.

#> Layer name: nc
#> Metadata:
#>   DBF_DATE_LAST_UPDATE=2016-10-26
#> Geometry: Polygon
#> Feature Count: 100
#> Extent: (-84.323853, 33.881992) - (-75.456978, 36.589649)
#> Layer SRS WKT:
#> GEOGCRS["NAD27",
#>     DATUM["North American Datum 1927",
#>         ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["latitude",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["longitude",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     ID["EPSG",4267]]
#> Data axis to CRS axis mapping: 2,1
#> AREA: Real (24.15)
#> PERIMETER: Real (24.15)
#> CNTY_: Real (24.15)
#> CNTY_ID: Real (24.15)
#> NAME: String (80.0)
#> FIPS: String (80.0)
#> FIPSNO: Real (24.15)
#> CRESS_ID: Integer (9.0)
#> BIR74: Real (24.15)
#> SID74: Real (24.15)
#> NWBIR74: Real (24.15)
#> BIR79: Real (24.15)
#> SID79: Real (24.15)
#> NWBIR79: Real (24.15)

without directly calling system or ogrinfo.


texts: Framework/software for dynamic WPS creation

texts:
We have some specific EO data processing facility and are about to provide an access to one using WPS-based interface. The are plenty of solutions for the task. However, it is hard to find one that meets our requirements:

Ability to deploy services (or processes) automatically (without UI)
Ability to deploy service without rebuilding from sources or restarting the container. I.e. deployment of one service/process shouldn't affect others.
Open/free license

There also some Java frameworks that allow to implement WPS service right in the app, but they are poorly documented.

Could anyone recommend something that suits our needs?


texts: Editing EGDB data with Portal/Collector: Do error messages from the DB pop up in the app?

texts:
Oracle 18c; Portal/ArcGIS Server 10.7.1
I'm considering setting up Portal  so that I can edit my on-premise EGDB data in the field using ArcGIS Collector.
As part of this implementation, I would have triggers in the Oracle database that would synhronously integrate the EGDB data with an external system (via Oracle's UTL_HTTP package).
If a record were edited, but failed to sync, the Oracle trigger would return an error message.
If I were editing in ArcMap, then a message would pop up in a dialog and I wouldn't be able to save the record.
But I'm not editing in ArcMap, I'm going to be editing in Collector, so my question
is:
Would an Oracle database error pop up in a message box in Collector? And would the user be able to save the record?


texts: within polygon: why SpatialIndex slow down query

texts:
I create such simple demo table:

CREATE TABLE IF NOT EXISTS tbl1(id INTERNAL PRIMARY KEY UNIQUE, name TEXT NOT NULL);
SELECT AddGeometryColumn('tbl1', 'geom', 4326, 'POINT', 'XY');
SELECT CreateSpatialIndex('tbl1', 'geom');

after that I create 180 * 360 * 60 points,
and do query:

spatialite> select count(*) from tbl1 where ST_Within(tbl1.geom,
PolyFromText('POLYGON((-179 53, 179 53, 179 50, -179 50, -179 53))', 4326)) == 1;
32399
CPU Time: user 1.552150 sys 0.063322

Too slow, then I want to use index:

spatialite> select count(*) from tbl1 where ST_Within(tbl1.geom,
PolyFromText('POLYGON((-179 53, 179 53, 179 50, -179 50, -179 53))', 4326)) == 1 AND
tbl1.ROWID IN (select ROWID from SpatialIndex where f_table_name = 'tbl1' AND f_geometry_column='geom' AND search_frame = tbl1.geom);
32399
CPU Time: user 3.325508 sys 0.066661

so instead of speedup, I have ~2x slow down, why?

Update

for the first one query explain query plan tells 0|0|0|SCAN TABLE tbl1

for the second one query:

0|0|0|SCAN TABLE tbl1
0|0|0|EXECUTE CORRELATED LIST SUBQUERY 1
1|0|0|SCAN TABLE SpatialIndex VIRTUAL TABLE INDEX 1:


texts:
I have question: is it possible to set up arcpy and pyqgis on the same IDE (I use Pycharm).
I have some problem with running my pyqgis scripts from pycharm. When I first configure Pycharm to recognize pyqgis (I set pythonpaths) everything works well, but after program restart I got this: ImportError: No module named qgis.core Is there any solution to this?


texts: Extraction not empty data featureclass from a workspace by arcpy!

texts:
What I want from python script,  searching a workspace containing Geodatabases and feature classes(Iterate Feature Classes), then clip those by interested Area (polygon). the resulted clipped files most not empty. I tried to convert Model builder to Python script and added "If" condition to transfer not empty features to final geodatabase:

import arcpy
from arcpy import env
#script arguments
env.workspace=arcpy.GetParameterASText(0)
clip_Features = arcpy.GetParameterASText(1)
FeatureClass_clip = arcpy.GetParameterASText(2)

#Local variables:
Recursive = "true"
Name = "FeatureClass"
Feature_Class = "\\FeatureClass"
FeatureClass_Clip__2_= FeatureClass_clip

FeatureClass = arcpy.ListFeatueClasses("*")
#Process : Clip
arcpy.Clip_analysis(FeatureClass, Clip_Features, FeatureClass_Clip, "")
#Process: Delet Features
for fc in FeatureClass_Clip:
    Count1=str(arcpy.Getcount_managment(fc))
    if count1 <> "0":
       arcpy.FeatureClassToGeodatabase_conversion(FeatureClass_Clip, "")

I imported the script to toolbox and defined parameters, but script
not ran because of clip raise error.


texts: Printing map portions one by one in QGIS?

texts:
I am using QGIS 2.14.2 and I have a map divided into different regions. Each piece have a number assigned (eg. 27). I have a page configured for printing the region (territory) along with 4 pieces of attribute tables filtered by a field related to the number (eg. 27) corresponding to the map to display information related to the area (eg. a territory street list with the 'territory' field equal to 27).

How to set up the easy way to print the other territories; for example, print the territory number 26, the 25, the 15, or another one? ...

... And allow the table data (eg. a street list of the territory) also present the information related to that territory?

I need to configure one by one or there is a better way?


texts: GeoServer: `org.geotools.process.raster.AffineProcess` could not be instantiated, how to fix?

texts:
I have installed GeoServer from OpenGeo Suite 4.1 on a fresh Unbuntu 12.04 following Boundless instructions. Now I wanted to test geoscript-groovy, and installed groovy, geoscript-groovy-1.3.1, geoserver-2.5.2-wps-plugin and geoserver-2.5-SNAPSHOT-groovy-plugin (community, 10-15-2014), with:

groovy -v
Groovy Version: 1.8.6 JVM: 1.6.0_32 Vendor: Sun Microsystems Inc. OS: Linux

When I restart tomcat 6, I got the following warning about Errors that a class called org.geotools.process.raster.AffineProcess cannot be instantiated.

Does anyone know what's the cause of this error, and how to fix it?

  catalina.2014-10-16.log:WARNING: Can't load a service for category
  "RasterProcess". Cause is "ServiceConfigurationError:
  org.geotools.process.raster.RasterProcess: Provider
  org.geotools.process.raster.AffineProcess could not be instantiated".

  catalina.2014-10-16.log:java.util.ServiceConfigurationError:
  org.geotools.process.raster.RasterProcess: Provider
  org.geotools.process.raster.AffineProcess could not be instantiated

  catalina.2014-10-16.log:Caused by:
  java.lang.ExceptionInInitializerError catalina.2014-10-16.log:WARNING:
  Can't load a service for category "RasterProcess". Cause is
  "ServiceConfigurationError: org.geotools.process.raster.RasterProcess:
  Provider org.geotools.process.raster.AffineProcess could not be
  instantiated".

  catalina.2014-10-16.log:java.util.ServiceConfigurationError:
  org.geotools.process.raster.RasterProcess: Provider
  org.geotools.process.raster.AffineProcess could not be instantiated
  catalina.2014-10-16.log:Caused by:
  java.lang.ExceptionInInitializerError


texts: How to have a record id without geometry?

texts:
I have a buildings shapefile layer and every building has a building ID linked to its feature id ($id). Some new buildings have been built and old buildings were demolished. I'd like to not recycle those old ID's, I'd like to sequentially add new numbers to the new buildings, but I would like to delete the geometry of the old buildings.

Is it possible to have a record in my database that does not contain geometry?

Is there a better way to do this?

I mainly use QGIS but I have the basic version of ArcGIS Desktop and ArcGIS Pro at my disposal if needed.

I'm somewhat familiar with PostgreSQL and geodatabases if those are recommended over shapefiles.


texts: Printing a map using GeoExt, Openlayers and GeoServer

texts:
I am trying to print an Openlayers map using GeoExt and Geoserver. Here is the code

printProvider = new GeoExt.data.PrintProvider({
        url: printServiceURL,
        method: "POST", // "POST" recommended for production use
        capabilities: printCapabilities, // provide url instead for lazy loading
        listeners: {
            "print": function () { printDialog.close(); },
            "printexception": function(printProvider, response) { alert('Error ' +    response.statusText); }
        }
   });

I am getting the desired output when I am using GET method but it does not work using POST method. What could be the issue ?


texts: Separating spectral similar training sets to improve accuracy of classified Landsat scenes in ArcMap?

texts:
I am classifying Landsat TM scenes of the entire Chesapeake Bay as part of my Master's research in ArcMap 10.3. For these scenes, I am classifying the following land cover types: water, forest, agricultural, wetland, beach, and urban. Training sets have been drawn throughout the scene with water having the highest count of pixels (2731803). Beach has the least number of pixels (2944) but is my most important cover type and there is no more available beach to delineate training sets on.

As one would expect, my urban and beach classification categories are extremely spectrally similar and make classifying these land cover types particularly difficult. Looking at the scatterplot comparing both categories, the two classes overlap heavily. An example can be seen below; red is beach and purple is urban:

Typically, categories that overlap in the scatterplot are merged to create a single training set but that is simply not possible to do for my research. I have tried adding more training sets for both cover types, merged and deleted repetitive training sets, and tested different parameters within the Maximum Likelihood classification scheme all to no avail.

An example of a misclassified island can be seen below. These occur throughout my entire scene. The red pixels in the image are misclassified beach pixels.

How best can I go about separating all of my land cover classes and reducing the frequency of misclassified pixels in my scenes?

I have been unable to find papers and posts regarding this issues.


texts:
Anyone know of a free OpenSource .NET library i could use to convert KML to WKT?

I know i can use openlayers and there are many JavaScript possibilities, but i need to do this on the back-end.


texts:
I have a polygon shapefile with 800 records. What I want to do is to create an MXD document and a map in jpg for every polygon record from the shp.

Is this posible to be done with data driven pages or map book in ArcGis or by using any other tool?


texts: Seeking tool to calculate "one second of longitude equals how many meters" based on anywhere in world?

texts: I have 4 text boxes that are data driven based on "1 Sec Long Meters S". The problem is this field is not filled out. I'm wondering if there is a tool that will calculate how many meters one second of longitude is based on my location. I will also need to do this for latitude.


texts: Capture any embedded map larger than screen

texts:
I use the various map portals on the daily basis, like below:

and the biggest problem is, that my capture is restricted to the screen, being precisely to the area, where the map is embedded, as you can see above.
I found several tools, which are able to capture the maps larger than screen, but unfortunately they work only for Google Maps.
https://superuser.com/questions/145642/how-to-take-a-screenshot-of-a-something-bigger-than-the-screen
https://www.bpwebs.com/print-a-large-google-map-larger-than-the-screen/
https://stackoverflow.com/questions/22210034/making-screenshot-of-google-map-bigger-than-my-device-screen
https://webapps.stackexchange.com/questions/5734/can-i-extract-a-google-maps-mashup-bitmap-that-is-larger-than-my-screen
I need something which could be able to grab the map content from any existing interactive map (i.e. LEaflet, Openlayers, and so on) with a resolution larger than the screen.
Is it possible at all?


texts: Display the size of the dot based on the value

texts:
I import a CSV File and for each coordinate, I get a point in the map.
Can I set a column for the value. Then the size should based on the value have a special size. Or based on all values, it should display from max to min value.


texts: gdal.Warp() returns different results on 2.1.0 and 2.3.2

texts:
I upgraded GDAL version from 2.1.0 to 2.3.2. Both versions I built from source with Python bindings. PROJ version - 5.2.0.

Expected behavior

gdal.Warp() returns the same result as in 2.1.0

Actual behavior

gdal.Warp() returns result with some changes in raster array.

Steps to reproduce the problem.

The archive contains source.tiff.
rasters.zip

I have the code:

import hashlib
from osgeo import gdal

source_filename = "source.tiff"
gdal_version = gdal.__version__
warp_options = '-of GTiff -r cubic -te_srs "+proj=longlat +datum=WGS84 +no_defs " -te 49.21875 -14.5197800463 49.306640625 -14.4346802153 -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" -ts 256 256 -dstalpha -srcnodata 0 -wo UNIFIED_SRC_NODATA=YES'

def create_output_filename():
    return "processed_{v}.tiff".format(v=gdal_version)

def warp():
    output_filename = create_output_filename()
    gdal.Warp(output_filename, source_filename, options=warp_options)

warp()

I run it on GDAL 2.1.0 and GDAL 2.3.2.
The results are two .tiff files - processed_2.1.0.tiff and processed_2.3.2.tiff
The archive contains results, too.

How to check diff: 1) gdalinfo:

gdalinfo processed_2.1.0.tiff -stats

...
Band 1 Block=256x8 Type=UInt16, ColorInterp=Gray
  Minimum=5841.000, Maximum=26698.000, Mean=7100.601, StdDev=1761.655
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=26698
    STATISTICS_MEAN=7100.6005859375
    STATISTICS_MINIMUM=5841
    STATISTICS_STDDEV=1761.6546349845
Band 2 Block=256x8 Type=UInt16, ColorInterp=Alpha
  Minimum=255.000, Maximum=255.000, Mean=255.000, StdDev=0.000
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=255
    STATISTICS_MINIMUM=255
    STATISTICS_STDDEV=0

gdalinfo processed_2.3.2.tiff -stats

...
Band 1 Block=256x8 Type=UInt16, ColorInterp=Gray
  Minimum=5841.000, Maximum=26699.000, Mean=7100.600, StdDev=1761.724
  Metadata:
    STATISTICS_MAXIMUM=26699
    STATISTICS_MEAN=7100.5995788574
    STATISTICS_MINIMUM=5841
    STATISTICS_STDDEV=1761.7243529314
Band 2 Block=256x8 Type=UInt16, ColorInterp=Alpha
  Minimum=65535.000, Maximum=65535.000, Mean=65535.000, StdDev=0.000
  Metadata:
    STATISTICS_MAXIMUM=65535
    STATISTICS_MEAN=65535
    STATISTICS_MINIMUM=65535
    STATISTICS_STDDEV=0

2) check raster arrays:

2.1.0:

2.3.2

Operating system

Ubuntu 16.04 64 bit

GDAL version

2.1.0 and 2.3.2, built from source

Any ideas?


texts:
I have been working with the atlas generation option in QGIS for some time and it works great! What I cannot figure out is how to incorporate the overview function into it.

Here is an example of what I want to get. I want to generate an atlas of 10 different UTM grids. One map per grid. On each map I want a zoom-in overview map to the shapefile labeled with numbers on the map.

However when I move on to the next UTM grid in my atlas generator, the overview is still on the old location. Here I would like an overview map showing number 17 to 19, not 9 and 11, since they are no longer in the scope of this map.

If I make the overview map to be selected by the atlas I can only have it selected by the same UTM grid I have for my main map, which is off course not very usefull. Any thoughts on how to get this done? Feels like I should be able to do some conditional formatting on the atlas generation level.


texts:
I'm trying select features from one shapefile that contains features from another shapefile. However when I run the command I get an error " no such column geometry" Both shapefiles are working normally.
Any ideas what could be the issue?
Here is my command:

ogr2ogr -f "ESRI Shapefile" selection_polygon.shp shapefile_polygon.shp -dialect sqlite -sql "SELECT * FROM shapefile_polygon polygon, 'shapefile_point.shp'.shapefile_point point WHERE ST_Contains(polygon.geometry, point.geometry)


texts: Plot a Raster in R with white space for NA values

texts:
I want to produce a .png of a raster with background values plotted as white space.

Below you can see the current state of the plot.

As you can tell, the pixels which should contain NA values, in fact register as a single value of 94.  The issue is that running a simple replace via r[ r== 94] <- NA would also conceivably replace points within the raster.  Also it's a 21GB raster, so this method would be computationally onerous.

I've produced four other .pngs from the same raster stack without this issue, as png(raster.png); plot(raster); dev.off() automatically produced white space for background values.  Do you have any other suggestions for dealing with this, outside of a crop or something in photoshop?


texts: Ignoring group layers when iterating through layers

texts:
I have this nice code that goes through each dataframe in the current MXD, and through each layer in the dataframe, and generates a list containing the name of the dataframe, the name of the layer, and the name of the coordinate system of that layer, divided by commas. Something like this:

Focus Area,SubArea,South_Pole_Lambert_Azimuthal_Equal_Area
Focus Area,SubArea_2,South_Pole_Lambert_Azimuthal_Equal_Area
and so on.

The code is this:

import arcpy
mxd=arcpy.mapping.MapDocument("CURRENT")
dataframes=arcpy.mapping.ListDataFrames(mxd)
for dataframe in dataframes:
    for lyr in arcpy.mapping.ListLayers(mxd, "", dataframe):
        print(dataframe.name+","+lyr.name+","+arcpy.Describe(lyr).spatialReference.name)

The code works well if there are no Group Layers in the dataframes. If there is a group layer containing layers, the code will stop and throw this message:

Runtime error
Traceback (most recent call last):
  File "<string>", line 6, in <module>
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\__init__.py", line 1234, in Describe
    return gp.describe(value)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 374, in describe
    self._gp.Describe(*gp_fixargs(args, True)))
IOError: "Climatology and environment (L)" does not exist

So, "Climatology and environment (L)" is a group layer. In my understanding, the code is trying to find the "layer" and get its coordinate system, but there is no actual layer, so it just breaks. Do you think of any workaround to this? A way to ignore the group layers maybe?


texts: IGPFunction2.Execute() never triggers when using IEnumLayer.Next() when run in background

texts:
I got a strange issue in GP tool I don't really understand:

In my IGPFunction2.UpdateMessages() method, I use the following code to traverse all layers:

IEnumLayer enumLayer = gpUtils.GetMapLayers();
enumLayer.Reset();
ILayer layer = enumLayer.Next();
while (layer != null)
{
    <code here>
    layer = enumLayer.Next();
}

The purpose is to look for layers with duplicate names.

This works fine when I don't run the tool in the background. In background mode however, the Execute()
method is never triggered.

This seems to happen after enumLayer.Next() is used. When removing the enumLayer.Next() and the code below, it works fine.

There are no error messages or exceptions occuring, and Arc claims that the tool executed succesfully.
When calling IGPMessages.AddMessage() as the first action in Exceute(), nothing is yielded.

Any suggestions is appreciated..


texts: Does anyone of you know how and where I can download postal area districts of Germany in shapefiles (download free of charge if possible).


texts: DEM/Raster Peak and Valley identification and/or elimination (non-generalization)

texts:
Is there is a tool or process that can identify and/or eliminate obvious peaks and valleys in a DEM?

I have looked into smoothing at  What raster smoothing/generalization tools are available? but I am interested in keeping the original sharpness of the DEM.

The DEM is created using the topo to raster from a point cloud, so I seek suggestions on ways to identify the peaks and valleys from the pre-DEM point cloud.

I've looked at  How to calculate peaks in a DEM too, but I am interested to see other ideas, and ideas from an arcgis point of view.


texts: Using multiple St_Dwithin clause result in a huge performance loss

texts:
In PostgreSQL using PostGIS, I want to find every point which are close enough of some other points but which also are far enough from another points series.

So here is my query :

SELECT pod.id, st_AsText(point)
FROM
    pod,
    (SELECT poi.geom FROM poi WHERE poi.sector_id IN ('39')) AS restric,
    (SELECT poi.geom FROM poi WHERE poi.sector_id IN ('25')) AS exclu
WHERE ST_DWithin(restric.geom, pod.point, 10)
AND NOT ST_DWithin(exclu.geom, pod.point, 5)
GROUP BY pod.id

At so point my query will return me all the point I am interested in but it takes multiple minutes to answer.
Whereas when I request only one of the two clauses of my query, PostgreSQL answer in no time (less than a second).

Is there a reason for that ? And is there a way to keep the execution time reasonably low ?

Here is the EXPLAIN (ANALYZE, BUFFERS) result :

QUERY PLAN
Group  (cost=5061.73..5830.65 rows=409 width=37) (actual time=7726.053..7820.974 rows=596 loops=1)
  Group Key: pod.id
  Buffers: shared hit=1892516, temp read=1799 written=1799
  ->  Sort  (cost=5061.73..5062.75 rows=409 width=37) (actual time=7726.032..7794.860 rows=328949 loops=1)
        Sort Key: pod.id
        Sort Method: external merge  Disk: 14376kB
        Buffers: shared hit=1892516, temp read=1799 written=1799
        ->  Nested Loop  (cost=2.72..5043.99 rows=409 width=37) (actual time=0.773..7616.931 rows=328949 loops=1)
              Join Filter: ((NOT (poi.geom && st_expand(pod.point, '0.130650640188137'::double precision))) OR (NOT (pod.point && st_expand(poi.geom, '0.130650640188137'::double precision))) OR (NOT _st_dwithin(poi.geom, pod.point, '0.130650640188137'::double precision)))
              Rows Removed by Join Filter: 43
              Buffers: shared hit=1892516
              ->  Nested Loop  (cost=2.72..11.15 rows=1 width=37) (actual time=0.585..3.602 rows=596 loops=1)
                    Buffers: shared hit=216
                    ->  Index Scan using "POS_pkey" on pos  (cost=0.42..4.44 rows=1 width=32) (actual time=0.026..0.028 rows=1 loops=1)
                          Index Cond: ((id)::text = '6'::text)
                          Buffers: shared hit=4
                    ->  Bitmap Heap Scan on pod  (cost=2.30..6.70 rows=1 width=37) (actual time=0.552..3.399 rows=596 loops=1)
                          Recheck Cond: (point && st_expand(pos.geom, '0.653253200940685'::double precision))
                           Filter: ((pos.geom && st_expand(point, '0.653253200940685'::double precision)) AND _st_dwithin(pos.geom, point, '0.653253200940685'::double precision))
                          Rows Removed by Filter: 268
                          Heap Blocks: exact=201
                          Buffers: shared hit=212
                          ->  Bitmap Index Scan on pod_geom_gist  (cost=0.00..2.29 rows=2 width=0) (actual time=0.404..0.404 rows=864 loops=1)
                                Index Cond: (point && st_expand(pos.geom, '0.653253200940685'::double precision))
                                Buffers: shared hit=11
              ->  Seq Scan on poi  (cost=0.00..4922.41 rows=409 width=32) (actual time=0.047..12.350 rows=552 loops=596)
                    Filter: ((sector_id)::text = '25'::text)
                    Rows Removed by Filter: 139241
                    Buffers: shared hit=1892300
Planning time: 0.650 ms
Execution time: 7823.665 ms


texts: Determining average values in single flowline using ArcGIS for Desktop?

texts:
I need to get the average (Velocity) values of 1000 points on my long flowline. I tried to make parallel lines on each gridpoint using the "copy Parallel" function. But there two problems. It takes too much time, and at the end I can't automatically construct points for all the perpendicular lines to the flowline. I need to do one by one.

How can I automatically construct points for all the lines? Or is there any other technique that I could apply (create a grid and later to create my flowline?)

Details

The flowline is 100km and has 1000 points. On each of them I need to have the width-average value. The width is 5km. In the center, where is the flowline, I have high Velocity values. Near to the sides the Velocity decreasing. I need to have the average value on each point from side to side. So, I create a "5km" line applied on each point perpendicular to the flowline. I found a way with a tool to construct points on each of 1000 small ("5km") lines. But, this takes too much time, and I need to have curves not only lines.


texts: Is it perfectly acceptable to produce a map at a scale of 1:24,000 from a database that was digitized from a 1:100,000 scale source map?


texts: ArcGIS ModelBuilder looping wrong tools?

texts:
I've created a toolchain in ArcGIS ModelBuilder (see in the picture)

But after Creating the "Flow Direction" and the "Iterate field Values" it starts again with the "DEM reconditioning", it's an endless loop. But why it is doing so and how can I stop that.

My target after creating the flow direction is to project the resulting layer to all UTM stripes like described in another question: Looping Projection in ArcGIS ModelBuilder?


texts: Symbolizer with external GIF image is not rendered pixel-perfect

texts:
We want to pinpoint some locations on a map with some external GIF symbolizers. The issue we have is visible in the attached image*: looks like somethings happens on a sub-pixel level, so that the resulting pluses look mutated :), most of them are somehow distorted/blurred/interpolated. The original external GIF is shown in the right lower corner of the image* (simple and sharp 5x5 pixel "+" symbol).

Thank you in advance for any clue or solution.

Regards,

skipr

Original image:

Magnified image:

Update: Same thing happens with well known symbolizers (e.g. MARK_CROSS):

...
final Mark mark = sb.createMark(StyleBuilder.MARK_CROSS, Color.red, Color.yellow, 1);
final Graphic graphics = sb.createGraphic(null, mark, null, 1, 10, 0);
final Style style2 = sb.createStyle(sb.createPointSymbolizer(graphics));
...


texts:
I have a shapefile with several polygons, and I want to convert it to a raster file with 30mx30m (using QGIS 3.2).

I am using "Raster > Conversion > Rasterize", and I am having trouble with the values I have to use.
So, for the "Input layer" I am obviously using the shapefile. In the "Field to use for a burn-in value", I am using a column which is filled with "1" values for all polygons.
For the "Output raster size units", I am selecting "Pixels", and both "Width" and "Height" is 500.000000 (the pixel size seems to get smaller each time this value is increased, and I don't fully understand why). "Output Extent" has the same extent as the shapefile.

This is how part of the maps looks like (using the parameters above):

And these are the properties I get in the file:

First question: I don't understand how the "Width" and "Height" values work; how do I set pixel size as 30m?

Second question: How do I get each dark pixel to be, for instance, "1", and the white pixel to be "0" (in a way that the raster is binary)? Or is this automatic?


texts: Number of intersected features from one layer per feature from second layer

texts: How can I get number of intersecting features from one layer per feature from second layer using any component from QGIS 2.8 (GRASS, SAGA, Python scripting...)? Result should be stored in new attribute column.


texts: Network Analysis Service Area Algorithm QGIS 3.16

texts:
I've been using ORS Tools to run my service area analysis for my research, because it was easier and quicker for me to get an isochrone, instead of the standard network analysis tools from QGIS, that needed me some extra steps to turn the road network into an isochrone.
I'm really curious on how it works. Are there any source I can use to study on the algorithm or the coding for this? I've tried looking for the algorithm behind ORS tools, but i didnt find any yet. Best I got was https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html for the, I'm assuming, coding behind service area analysis from the standard QGIS tools. I'm basically blind towards coding, Python, and stuff. I'm hoping to get a brief surface level explanation on how things work to get the end result.


texts: OpenLayers 3: How to check if vector source is ready?

texts:
ol.source.getState() doesn't seem to be reliable. When I call it on a vector source it returns ready, but the features are not available yet. Code looks like this:

var vectorSource = new ol.source.Vector({
  url: 'world.topo.json',
  format: new ol.format.TopoJSON()
});

// ... init map with vectorSource

console.log(vectorSource.getState()); // returns "ready"
console.log(vectorSource.getFeatureById("US")); // returns null

Any other way to see if a vector source is ready?


texts: Creating Multi-Part Polygons that touch internally using ArcGIS Desktop?

texts:
I have a set of parcels as illustrated in the image below. There are currently no multi-part features in dataset. What I would like to do is to merge (not the right word) all of the parcels with same address, say 421 for example. I want to have these three polygons as one multi-part polygon, that still has the same geometry, as in it doesn't dissolve the lines.

When I have tried this thus far, the lines dissolve resulting in only the boundary being preserved.

I feel like there is a super easy was to do this, but I can't find it. I have full access to ArcGIS Advanced and use it to manage the topology of this dataset.


texts: Substracting raster values from 2 rasters with different cell sizes

texts:
I am working with two raster files. The first is a DEM and the second an elevation-correction file. Both files display the same geographic region. My problem is that the cell sizes of both files aren't equal. The DEM contains 5 million pixels (values), whereas the second file contains only 30 which means the cell size is much broader. My wish is to subtract the second file from the DEM. I tried using Raster Calculator but the result is not satisfying since this destroys the whole DEM.

My approach is the following:
I have to use an expression within the Raster Calculator that tells the software to subtract every single DEM-cell from the corresponding value of my correction file. But unfortunately I haven't figured out how to do this...
Does anyone of you have an idea how to do this?

I am using ArcGIS 10.5 and QGIS 3.4.


texts: Converting GPS points in WGS84 to local Still Water Mine Grid in ArcGIS Desktop?

texts:
How to convert GPS points in WGS84 to local Still Water Mine Grid?

I got GPS points in WGS84, and just have info: Datum Still Water 0.0  to NAD_27_Montana South: 1903090, 506795 , rotation = -19.9988.

Any ideas about this?


texts:
I am looking for a way in QGIS (1.8.0) to visualise data, which is related to an area (polygon on a shape file). There are no x-y-coordinates with the data, just the info of the belonging area, means at an arbitrary position inside the area.

Example:
Samples x1,x2 and x3 are located in area a1, y1 to y5 in area a2.

As the samples have certain properties on which I intend to do queries on, distinct locations seems somehow necessary (e.g. sample x1 and x3 have a certain property and therefore should be shown/highlighted ).

I initially thought of assigning each sample some random coordinates from inside the areas, but what a hassle..

So, how to approach that??


texts:
I'm trying to show an array of selected field values in HTML in the map composer.

So far, this works but only for a specific value1 :

[% array_to_string( array( attribute(get_feature('layer', 'attribute1', value1), 'attribute2')))%]

What I want is to collect all the attribute from multiple selected features. I tried using is_selected() without success so far.


texts: I'm just doing a basic topograph and the interpolate line tool does not work. The DEM is selected and I even made a hillshade to see if that would solve the problem. After I finish the line and double click the line disappears. There is no error message, just a sound as if I was clicking on the wrong window but there are no other windows open. I've restarted GIS, I've made a new map, I've changed where the topo data goes, I've deactivated and reactivated 3D analyst and I've changed the coordinate system to a new one and back to the original and it still does not work.


texts: How to summarize overlapping polygon values within a single layer?

texts:
At ArcGIS Desktop 10.5, I have a polygon file geodatabase feature class.  Many of its polygons are spatially coincident; that is, if I click on a single polygon in Arcmap to select it, the resulting selection will have multiple records, all sharing the same exact polygon outline.  Here's a screenshot example where one polygon is selected, resulting in 74 selected records.  In my mind, they represent 74 polygons stacked upon each other.

The attribute table contains a long integer field, noise_value.  I need to sum this field for each input polygon, which, as I mentioned above, may be comprised of multiple spatially coincident polygons, and output the result to a new, flattened, polygon feature class containing the summed value.

In my example, the 74 overlapping polygons would be reduced to a single output polygon with a summed noise_value.  Of course, this summation would have to be repeated for all input polygons.

I've tried these approaches, and all have failed:

1. Flatten the input feature class while applying the field map merge rule "SUM" for the noise_value attribute.

I tried this approach with the Feature Class to Feature Class tool as well as the Merge tool.  In each case the output was identical to the input; no flattening or summing occurred.

2. Create a flattened output feature class and then spatially join it to the original feature class, summing by noise_value with field map merge rule SUM.

I tried this approach with the Feature to Polygon tool, and also with the Union tool to create the flattened output feature class.  Both cases failed, with a "Topoengine error".  Research indicates that this type of error can be traced to inadequate hardware resources, but I'm using a computer with an i7 processor, 16GB ram and a 1 TB SSD.

3. Use the Dissolve tool (@BERA's suggestion) with noise_value as the statistics field and Statistics Type = SUM.

Since what I want to do is flatten (dissolve) the geometry, I left the Dissolve Field blank.  The result did not work for me because it created a single output polygon with one overall summed noise_value.

I've run out of ideas, and am looking for advice.


texts: Selecting line inside polygon and copy attribute into them

texts:
I'm using QGIS 3.6.
I've got around 300 polygons with "name" field.
I would like to select lines inside those polygons and copy what's inside of "name" field into similar field of every line that is selected. Is there any smart method to do this and save hours of work using "Select by location" and copying name by hand?


texts: How does one calculate an angle in a geographic coordinate system?

texts:
Given three points on a line in a geographic coordinate system, how does one calculate the angle formed by these three points?

Please, give a general solution but also include system specific examples if you like.


texts: Is it possible to add an export button to Qgis2Web web interface after filtering certain layers?

texts:
I created a webmap with pointdata using Qgis2Web.
There are a few filters available, and I would like to give users the opportunity to export the filtered points to for example a CSV or XLSX.
Is that possible?


texts:
I am new in GIS programming, I have to created a geoportail  with the Configurable Map Viewer (CMV). I want to add my KML layers but they don't appear in the map. Below is the part of my code

operationalLayers: [{
            type: 'kml',
            url: "http://localhost:6080/arcgis/rest/services/ban/MapServer/",
            title: 'banques',
            options: {
                id: 'libelle',
                opacity: 1.0,
                visible: true,
                outFields: ["*"],
                mode: 0
            },
            editorLayerInfos: {
                disableGeometryUpdate: false
            },
            legendLayerInfos: {
                exclude: false,
                layerInfo: {
                    title: 'My layer'
                }
            }
  }, {
            type: 'kml',
            url:    'http://localhost:6080/arcgis/rest/services/ban/MapServer',
            title: 'Region Rabat Sale',
            options: {
                id: 'key',
                opacity: 1.0,
                visible: true,
                mode: 0
            }
  },

   ],


texts:
Interpreting 'GIS' in the geographical information science sense, I have a historical question about the 1936 OSGB datum following the 1936-62 retriangulation of Great Britain - hopefully not too off-topic..

Does anyone know of the reason for the 49°N 2°W "true origin" of the datum?  Ordnance Survey website's closest comment on this appears to be:

  The True origin latitude and longitude coordinates of the National
  Grid are 49 degrees north: 2 degrees west. The False origin latitude
  and longitude or coordinates are 49 degrees 45 minutes and 58 seconds
  north : 7 degrees 33 minutes 23 seconds west.

  The False origin which lies slightly southwest of the Isles of Scilly
  was devised to ensure that all National Grid coordinates were positive

So the "false origin" at 0 Easting 0 Northing (approx 7.56°W 49.77°N) is chosen so the Scilly Isles and western point of Northern Ireland would remain in positive coordinates - an admirably egalitarian nationalist justification. But 2°W seems a very England-centric longitude, approximately bisecting the 1.7°E and 5.75°W extremes of the nation.  And the only redeeming feature of 49°N seems to be that the 49°N 2°W "true origin" is marginally closer to a UK dependency (Jersey) than the French mainland.  So why this parallel/meridian intersection?


texts: ArcMap annotation "disappears" when moved to the edge of the projection in map frame

texts:
I have feature-linked annotations in a file GDB for a world map in Plate Carree projection.  When I edit annotations and move an annotation beyond the upper extent of the projection in the map frame (accidentally), the annotation "disappears."

It is still showing as "Placed" but now the fields for Shape Length and Shape Area are zero.  When I select the feature and Zoom to Selected, I am taken to a random place on the map.  I have tried panning or changing projection, but the issue seems to be that annotation feature is now infinitely small.

Unfortunately, edits were saved before this problem was realized.  Is there any way to fix this?

I am using ArcGIS Desktop v10.7.1


texts: Opening map package does not display Table of Contents?

texts:
When map packages (.mpk) that I have created with Arc 10.5.1 are subsequently opened, the map displays but the Table of Contents is blank:

But if I then click on any of the List by... buttons (other than List by Drawing Order), the Table of Contents now appears:

But List by Drawing Order will never display the Table of Contents:

This problem occurs with .mpks created with either the Package Map tool or the File > Share As > Map Package option.

How can I create .mpks where the Table of Contents automatically displays upon opening?


texts: I need to calculate the average daily temperature for each month. As an input I have 5 .asc grids for 6, 9, 12, 15, and 18  GMT for each month of the year. I was wondering if there is a predefined tool to calculate the average using polynomial fitting. There is one GRASS module r.fitpoly, however it is not available in the standard version. The information about the GRASS module I had found here: http://www.sciencedirect.com/science/article/pii/S1364815205001593. Could you please give me a hint what else I can use for such procedure?


texts: Cannot get WFS features from GeoServer using OpenLayers 3

texts:
I read a book about OpenLayers 3 and I guess its about an older OpenLayers version because it mentions a ol.source.ServerVector vector source for getting WFS layers from a server.

In the latest OpenLayers 3.9.0 that I use, when I search the API the only vector sources I can find are vector, cluster, tileVector and imageVector.

So my first question is, is ol.source.ServerVector deprecated? Should I use plain ol.source.Vector instead to fetch WFS features from my GeoServer service?

When I search "WFS" on the OL examples I get only one example that uses ol.source.Vector to get WFS from GeoServer.

My second question is about the settings.

I don't get what this does url: function(extent, resolution, projection) {, how the URL is implemented and what the strategy exactly does. Can somebody explain or provide links to tutorials?


texts: Augmented virtual reality with GIS capabilities

texts:
I want to build an AVR App for any platform or OS with GIS capabilities. Do you know a software for this purpose; I prefer something with less coding.

I found this source http://socialcompare.com/en/comparison/augmented-reality-sdks but it contains too many softwares and I don't know which is for what I want.
An example of what I would like to do you can find in the following website http://moptil.com/

This App should be built by using the camera of a mobile phone or a Tablet. Using this App one can see the real world and some spatial data and general information. My problem is that I don't know what software I should use.


texts:
I am running running Mac OS X (10.7.5), PostgreSQL (9.2) and PostGIS (2.0.2) (installed via Postgres.app).

The PostGIS extension seems to work perfectly with vector and raster layers. Once I enable PostGIS for the specific database I am able to create vector and rasters from the SQL Editor of phpAdmin3 (like with ST_MakeEmptyRaster).

The problems begin when I try to load a raster with raster2pgsql. I have no great experience with Postgres (or DBMS in general) so what I did might appear stupid...

I tried to launch it via psql following this steps:

1) Open psql terminal;

2) typed \c mydatabase;

3) Got as prompt mydatabase=# (that I guess means I am logged with owner rights...;

4) Now I tried raster2pgsql -Gbut of course is no SQL so I got ERROR:  syntax error at or near "raster2pgsql" LINE 1: raster2pgsql -G;

Then I went back to my Mac OS X terminal shell and into the the folder where my raster2pgsql file is located (/usr/local/pgsql/bin) and typed raster2pgsql -s 4236 -I -C -M *.tif -F -t 100x100 /my/path/to/the/file/myfile.tif myrastertable | psql -U myuser -d mydb -h localhost -p 5432 but I got -bash: raster2pgsql: command not found.

What I am doing wrong?


texts: Can't load a shapefile with Python Console in QGIS?

texts:
I'm trying to load a shapefile into QGIS using iface.addVectorLayer, but whenever I do this, I get an error message "Layer is not valid: etc ..."

However, if I try to load the same shapefile use the menu dialog "Add Vector Layer", it loads fine. I have loaded raster files, so I know it's not some funny with the path name.

I'm sure I'm missing something obvious here, but what?

This line works:-

layer1 = iface.addRasterLayer("H:\Mapping Stuff\NE1_50M_SR_W\NE1_50M_SR_W\NE1_50M_SR_W.tif", "basemap")

This doesn't:-

layer2 = iface.addVectorLayer("H:\Mapping Stuff\ne_50m_urban_areas\ne_50m_urban_areas.shp", "urban", "ogr")


texts:
I'm an archaeologist and have been using photogrammetry to record rock art sites in northern Australia. That is working very well (we're using Agisoft Photoscan, it's fantastic).

My question relates to taking 3D models into ArcGIS 10.3 for further work. The end point I am aiming for is to be able to create vector datasets that represent the art within a geodatabase. This, I'm hoping, can be done by 'tracing' evidence of artwork to create new records within feature classes within a geodatabase, together with domains and so on to classify artwork. GIS 101, really. Getting to that point is the problem, as the imported models lack textures or meaningful symbology in ArcScene.

Can anyone recommend a general workflow for doing this, or good resources to help? ArcGIS helpfiles are a little vague around these issues. At present, and in brief, I have been:

Exporting the model as a collada (*dae) file from photoscan.
Downsampling/editing/ in meshlab
Importing as a multipatch into ArcScene

That works fine, mostly, but how would one apply texture files? Also, others I have spoken to suggest exporting/creating LAS pointclouds, which also works, but again I am stuck with no meaningful symbology at the end of the journey.

As you can probably tell, the 3D modelling thing is new to me.


texts: Aligning labels with line connecting feature to its label in QGIS?

texts:
I have labelled a points layer by defining a line between points and labels, because points were overlapping too much. Using additional columns in the points layer attributes table (LABEL_X and LABEL_Y), I can move the labels manually.

To achieve that, I have added a 'Geometry Generator' in the Style section of the layer's Properties dialog. To define the geometry generator, I have used the following expression where $x and $y are the features' coordinates and LABEL_X and LABEL_Y are the labels' coordinates:

make_line(make_point($x,$y),make_point("LABEL_X","LABEL_Y")).

My issue is that labels do not align properly with the line connecting them to features. There seems to be a default offset. How could I correct it? The images below show (1) the way labels display inside the map canvas and (2) the Geometry Generator definition.

Image 1 : Labels with connecting-to-features lines

Image 2 : Geometry Generator definition


texts: Would you consider online geocoding a breach of privacy?

texts:
Suppose I have a bunch of addresses of individuals participating in a certain study (most likely - health related, where privacy and ethical considerations are always important issues).

Nowadays, providers like Google or Yahoo offer decent results in terms of positional accuracy.

The North American Association of Central Cancer Registries (NAACCR) lists such options in their 'Geocoding Best Practices: Review of Eight Commonly Used Geocoding Systems' and 'A Geocoding Best Practices Guide' guides.

Cinnamon and Schuurman (2010) for example used BatchGeocode service as a part of their tool to investigate injuries in low resource setting.

Would you consider geocoding such addresses using online services, like Google Maps or OpenStreetMap a breach of privacy?

A possibly related question is Geocoding USA addresses that cannot be sent over internet?

An article in Epidemiology (one of the leading, peer-review journals in the field) published short communication detailing instructions on how to geocode using Google Maps & Places APIs. Interestingly, not a word about security/privacy was mentioned.


texts: retrieving polygon shapes from postgis in geodjango

texts:
Can anyone explain how to extract polygon shapes that are stored in a PostGIS database and display it in the browser?

I am using geodjango framework in python. I would like a step by step guide on how to achieve this.


texts: Determining mathematical expression for Boolean logic with 5 raster layers in ArcGIS for Desktop?

texts:
I have 5 layers in ArcGIS 10.2 for Desktop, and i want to apply Boolean Logic. What will be the mathematical expression in "Raster Calculator".

I know that for to layers expression will be like ("layer1"==1)&("layer2"==1)
I am trying following expression for 5 layers all together ... but I am getting an error message.


texts: I received some files from a third party that are supposed to be geoPDF. Is there a way to open a file (programmatically with .NET or Python, with a text editor, etc.) and verify it was created as a geoPDF? Does a geoPDF have bytes or an internal structure we can check?


texts: Poor image quality of raster pyramid files in ArcGIS Desktop (no anti-aliasing)?

texts:
There are several GeoTIFF files (size ~7000x7000 px) which I try to open in ArcGIS Desktop 10.0. I have already created raster pyramids for them (.ovr files) using geoprocessing tool but their quality is rather bad (0.25 resize):

Each GeoTIFF file has many thin details - actualy it is CAD drawing saved to raster. I have tried to do manual resize for GeoTIFF in free image viewer and used resample + anti-aliasing. Result looks great (0.25 resize):

Standard resample methods (in ArcGIS 10.0) which I have used for pyramid creation didn't give me any good result. I also tried to view image without pyramids - same thing.

Then I have tried to create that pyramid file (.ovr) using GDAL library (GDALDataset.BuildOverview). But GDAL also is not capable of using anti-aliasing when generating OVR files.

Another way is to generate OVR file manually. I have tried to find any library to create OVR file. Nothing. Soon, I have realized that OVR actually is a multipage TIFF file. When I have created OVR file manually in image editor (as multipage TIFF), renamed and tried to use it in ArcGIS - I saw black screen with some distorsion. Some problem with color palette may be. I also have not found any OVR specification.

Does anybody know is there any possibility to generate such pyramid files using anti-aliasing to improve image quality? May be any links to OVR specification or libs?

EDIT1:

Here is part of raster which can be used as an example (try to build pyramids in ArcGIS).

EDIT2:

In order to completely answer to a comments concerning changing resampling method when building pyramids I have tried all methods in ArcGIS. Then I have opened each ovr file as multipage tiff and took one pyramid level and saved it to a separate file. Then I have resampled same source image in freeware viewer (IrfanViewer) to the same size (pyramid level). Source image has size 5001 x 3361 px. Pyramid level has size 313 x 210 px. So, image was reduced in 16 times. Here is my results:

Nearest resampling method (from pyramid file):

Bilinear resampling method (from pyramid file):

Cubic resampling method (from pyramid file):

Resample with anti-aliasing in IrfanViewer to the same size:

EDIT3:

I am providing one of my source georeferenced rasters and its resampled example to show the quality which I am trying to achieve in pyramids.


texts: raster2pgsql loads raster into a single column

texts:
I would like to load a raster (.asc) into my PostgreSQL Database.
In order to demonstrate the problem I chose the raster downloaded from this homepage: https://shop.swisstopo.admin.ch/de/products/height_models/dhm25200

Here some info about my PostGIS_Version:

psql -U postgres -c "SELECT PostGIS_Version()"

    postgis_version
---------------------------------
2.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

Let me give you some information about this raster with gdalinfo:

Driver: AAIGrid/Arc/Info ASCII Grid
Size is 1926, 1201
Coordinate System is `'
Origin = (479900.000000000000000,302100.000000000000000)
Pixel Size = (200.000000000000000,-200.000000000000000)
Corner Coordinates:
Upper Left  (  479900.000,  302100.000)
Lower Left  (  479900.000,   61900.000)
Upper Right (  865100.000,  302100.000)
Lower Right (  865100.000,   61900.000)
Center      (  672500.000,  182000.000)
Band 1 Block=1926x1 Type=Float32, ColorInterp=Undefined
  Min=193.000 Max=4445.019
  Minimum=193.000, Maximum=4445.019, Mean=1272.461, StdDev=800.094
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=4445.019
    STATISTICS_MEAN=1272.4613092273
    STATISTICS_MINIMUM=193
    STATISTICS_STDDEV=800.09392403263

In order to load this .asc file I use this command in the command line:

raster2pgsql  -I -C -M "C:\Users\user\Downloads\data\data\DHM200.asc" -F public.dhmtest3 | psql -U postgres -d GIS_DB

To get this response:

BEGIN
CREATE TABLE
INSERT 0 1
CREATE INDEX
ANALYZE
HINWEIS:  Adding SRID constraint
HINWEIS:  Adding scale-X constraint
HINWEIS:  Adding scale-Y constraint
HINWEIS:  Adding blocksize-X constraint
HINWEIS:  Adding blocksize-Y constraint
HINWEIS:  Adding alignment constraint
HINWEIS:  Adding number of bands constraint
HINWEIS:  Adding pixel type constraint
HINWEIS:  Adding nodata value constraint
HINWEIS:  Adding out-of-database constraint
HINWEIS:  Adding maximum extent constraint
 addrasterconstraints
----------------------
 t
(1 Zeile)

COMMIT
VACUUM

That means, my raster was loaded to the DB. But looking in the DB, the raster is composed of only one column, with all the single raster values coerced together. Which means, this is not a raster but something different.


texts: Doing spatial regression using big data when Spatial weight matrix is too large for R?

texts:
I am currently doing an analysis on a quite huge data set - 250,000 polygons. I have figured out the general Moran's I for a linear model, but when I was trying to do a spatial lag or spatial error model, the error occurs:

  Error: cannot allocate vector of size 445.6 Gb

I am using R for the analysis, and I tried to use a sparse matrix but only "LU" sparse method can succeed, and the result is not really good enough.

Any solutions for that, perhaps using spatial regression based on Hadoop? Our lab can spend some budget on that, but won't be large.


texts: Getting the x of the next entry on the atribute table

texts:
I've been struggling for the past hours with this situation, I need to get the X and Y coordinates of the next entry of the atribute table like the image shows, is there any tool or function in QGIS that I can use?

Or do I have to use Excel or make a mini program of my own to do it? I've tried the "$X_at" and some other functions and nothing.


texts: Zonal statistics with feature collection column (Google Earth Engine)

texts:
I have run into a very similar problem a few times while working in Google Earth Engine and never found a good solution.
I am working with raster data (continuous values) and a feature collection of ecosystems, and my goal is to summarize the raster values over a column within the feature collection. The ecosystems shapefile has a column that assigns each ecosystem to a biome (so there are 827 ecosystems and 14 biomes). I want to calculate the mean raster value over these biomes.
In ArcMap I would use zonal statistics as table and specify that I want to calculate the statistics over the biome column, but I have not found a good way to accomplish this in GEE. So far I have been using the reduceRegions function and specifying my feature collection as the collection, but this only lets me summarize the raster values over the 827 ecosystems.
Is there a function that would allow me to do this in GEE? Or would it be easier to just make the changes to the shapefile in Arc and then re-import it?


texts: How to create simple select script in ArcPy?

texts:
I am new to Python and want to create a script inside the ArcMap Python window. My current workspace is a series of rivers portrayed in the state of Ohio. I have a list of fish and the rivers (from my series of rivers) they are found in. Ideally, I would like the user to be able to enter the name of a fish and have the script highlight or select the rivers in the map.

For example, if I were a fisherman who was not from this area and wanted to catch a certain fish (such as a smallmouth bass) I could type the fish into the script as a raw input and have the map highlight the rivers where this fish could be caught.

I have an excel sheet with column titles as fish (i.e. largemouth bass, smallmouth bass, etc) and each row is a river name that they can be found in

How might I go about doing this?

Excel format:

Largemouth Bass (Col 1. Title)           Smallmouth Bass (Col 2. Title)
Ashtabula River                          Ashtabula River
Chagrin River                            Chagrin River
Cuyahoga River                           Conneaut Creek
Grand River                              Grand River
Huron River                              Huron River
Rocky River                              Rocky River
Sandusky River                           Vermillion River
Vermillion River


texts: Processing Model .model from Ubuntu to windows and "vice versa" Python Error

texts:
I have prepared a processing model in Qgis 2.10 under ubuntu 15.04, I tried to test it on a Windows(7) pc with Qgis 2.10 and I got a Python error message,
Then I tried the opposite way, I prepared the model on the Windows pc and transferred it to my ubuntu pc. Got the same error.
someone knows how to solve it?
i post the screenshot of the error


texts:
I have a script that iterates through a folder of jpgs, clips them to a grid, inserts them to an mxd, and exports them. My original problem was when I added the clipped raster from a standalone script to my mxd, it would automatically render it as a 1 band stretched image when a 3 band RGB is desired. In the python window I was able to use copy raster (management) with the "ColormapToRGB" parameter to create a layer that was displayed correctly. The problem now is when I run it from a standalone it doesn't automatically add the RGB to my TOC like it did from the python window. I instead have to add the RGB with another line like I did with my original clipped raster and I run into the same 1 band rendering issue. Anyone one know how I can get the RGB image to appear in my TOC with the correct display?

import arcpy

#Define workspaces
arcpy.env.workspace = "C:\\FolderPath"
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("C:\\FilePath.mxd")

#Define layer locations
df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
quads = "C:\\Geodatabase.gdb\\Quads_14"
qSelect = "C:\Geodatabase.gdb\qSelect"
outRast = "C:\\Geodatabase.gdb\\rastClip"

#Slice strings for SelectLayerByAttribute
qname = path[12:-4].replace("_", " ")
qnameUp = qname.upper()
sqlString = "NAME LIKE" + " '" + qnameUp + "'"

#Inserts each JPEG in folder and does all geoprocessing/map labelling/exporting
for jpg in arcpy.ListRasters("2*","JPG"):
    arcpy.MakeFeatureLayer_management(quads, qSelect, sqlString)
    desc = arcpy.Describe(qSelect)
    extent = str(desc.extent)
    arcpy.Clip_management(in_raster=jpg, rectangle=extent, out_raster=outRast, in_template_dataset=qSelect, nodata_value="256",clipping_geometry="ClippingGeometry", maintain_clipping_extent="NO_MAINTAIN_EXTENT")
    topo = arcpy.mapping.Layer(outRast)
    arcpy.mapping.AddLayer(df, topo, "TOP")
    arcpy.CopyRaster_management(outRast, out_rasterdataset="C:/FilePath", nodata_value="256", onebit_to_eightbit="NONE",colormap_to_RGB="ColormapToRGB", scale_pixel_value="NONE", RGB_to_Colormap="NONE")

arcpy.RefreshActiveView()
arcpy.RefreshTOC()
mxd.save()
del mxd, df


texts: Bulk insert into MS SQL Creating ObjectID's

texts:
I am running an ETL from a spatial view into publication database. The view has results in duplicate ObjectID's. I started working with the get_ids SP example from ESRI, but that seems to be built more for a single record insert and not a bulk insert. Here is my example...

DECLARE @id as integer
DECLARE @num_ids as integer
EXEC L.i16_get_ids 2, 1, @id output, @num_ids output;

INSERT INTO L.TAXPARCELWCONDOS ([OBJECTID]
           ,[PARCELNO]
           ,[SHAPE])
SELECT @id
      ,[PARCELNO] ,[shape]
  FROM [Parcel].[L].[view_TaxParcelwCondos];
END

How can I generate new ObjectID's for the insert.

Thanks,

Dave


texts: Reuse same style settings for chloropleth maps based on different datasets

texts:
I am new to QGIS. I have been using MapInfo for a couple of years, but I am considering to switch to QGIS. So far I like it, but there are some things that I haven't found out yet.

I want to make multiple graduated chloropleth maps based on a varying set of data but with the same graduated/rule based style.

For example, I have a polygon vector maps with postal code areas and would like to give the polygons a color based on the drivetime to a certain location A.

I can make this style as a graduated style or as a rule based style and save it.

Next, I need to make the exact same map for a different dataset, say drivetime to location B. I want to reuse the same settings (colors en cut-off points for the categories).

In mapinfo this was very easy:
Import an excel file with different datasets in each column-> make 1st map and save the thematic map settings -> Make a new map based on a new column but with the saved thematic map settings.

How do I do this in QGIS without
-having to remake the graduated style each time (saving a graduated style remembers which column was used and is gone when you change it)
or
-redefining all rules in the rule-based style?

Or is there a way to make a batch of maps that does work in QGIS 1.8? The tutorial on batching maps does not work anymore in 1.8 (http://www.qgis.nl/2012/07/28/series-kaarten-genereren/?lang=en )

Hope someone can help me.

THanks in advance,

niels


texts:
From techcrunch.

I see a lot of press coverage about Google Glass, yet I haven't seen details of GIS apps that use it.

What GIS applications use Google Glass?


texts: Extract just power features from OpenStreetMap

texts: Is there a way to only download power features from OpenStreetMap in some sort of geodata?


texts: PostGIS: Can´t preview geometry when in different schema than public - why?

texts:
So I'm trying to wrap my head around PostGIS and how to get it working correctly in a schema other than 'public'.
I followed this tutorial, and transfered PostGIS to my own schema, postgis_test. I used the PostGIS shapefile importer and imported my shapefile. Everything looks good (see screenshot of the geometry column), but the option to preview the geometry (the 'Eye'-icon) from the table is gone - why could that be?


texts:
While reconciling a child version with the parent version, in a dataset containing Geometric network, the following error was noticed.

Error during reconcile: "The version could not be reconciled. Invalid network element id : Field is not editable "

We work on ArcGIS 10.1 .

Tried to repair the Geometric network and do the reconcile , but to same effect.

What is the core issue behind the error and what is the workaround.

Experts pls advice.


texts: Is there a consumer-level handheld GPS that can use standard GIS data as its background layer?

texts:
I'm shopping around for a handheld GPS, and from what I'm seeing, none of the big vendors (Garmin, Magellan and Delorme) make a unit that can use, for example, a simple georeferenced TIFF or a Shapefile as a background layer on screen. They all want you to buy (or subscribe to) their own closed-off product which basically just repackages data that is often otherwise freely available (like USGS quads) and allows you to use it with their GPS units.

Is there a vendor that sells a device like a Garmin GPSMap or Delorme PN which can use freely-available GIS data standards like TIF, SHP, OSM data, and so on? I know Trimble's handhelds can do some of this, but they are much more expensive.


texts: Symbolising 2.5D polyline based on slope of each segment using QGIS?

texts:
I have a 2.5D polyline (each vertex has an elevation value) of a hiking route. I would like to change the line width depending on the slope of each segment in order to visualise how difficult each section of the route would be.

If necessary the segments could be split into individual features in a layer which could have a graduated symbology based on a calculated slope but ideally the route would remain intact as a single line.


texts: Translating latitude and longitude into "polar coordinates" around a point

texts:
Times Square in New York City has the coordinates (as per Google Maps):

40.7589, -73.9851

The Barrymore Theater has the coordinates (again, as per Google maps):

40.7601, -73.9862

Given that latitude and longitude are not rectilinear coordinates, how should I translate the coordinates for the Barrymore Theater into a radius from Times Square and an angle relative to due north? These would be like polar coordinates, with Times Square as the pole and due north as 0º.

I know that the Cartesian-to-Polar transformation is

radius = sqrt(x^2 + y^2)
angle  = tanh(y / x)

Is this a suitable approximation if I just plug in latitude and longitude for x and y? At what distance scale would it start to break down, and how badly?


texts: Converting lat-lon from degrees and decimal minutes to decimal degrees in PostGIS

texts:
Creating a geometry in Postgis is fairly simple:

ST_SetSRID(ST_MakePoint(longitude,latitude),4326);

These lat/lon should be in decimal degrees format like 33.02505 -96.70668

Is there a way/function to enter them in the degrees and decimal minutes format provided by NMEA $GPRMC sentences like 3301.5032 09642.4010?

Of course you can convert them in php or nodejs; but it would be easier (faster) to have it natively in Postgis.


texts: Using Configurable Map Viewer with Angular Application?

texts:
I have Angular JS application with WEB API RESET service
i tired to take sample from CMV Viewer and applied some feathers like search , print, legend  and included many layers from Our oracle DB

my case   :
I developed search screen in Angular JS application and want to display results on map

so how can using CMV viewer in my Angular application if someone use it before or any suggest ?


texts: Display coded value domains instead of descriptions in attribute window

texts:
Is it possible to display coded value domains instead of descriptions in attribute window in ArcMap?

I know this article shows how it can be done for the Attribute Table, but what about the Attribute Window?


texts: Leaflet- How to use time dimension plugin with image overlay

texts:
I am very new at Javascript.

I am trying to create animation with my monthly soil moisture output (.png format) on OpenStreetMap using the Leaflet package. I could use the time slider plugin to move images but I also want a play button that is available in the 'TimeDimension' plugin.

I have 150 images to display. I looked into the demo but I could not understand where to start as a beginner.

Is there any simplified example or instruction of using the plugin?


texts: Add X/Y coordinate format incompatible with viewshed tool coordinate format in QGIS

texts:
I am performing an r.viewshed analysis (id: grass7:r.viewshed) in QGIS 3.10.9 with Grass 7.8.3.
My question pertains to the 'coordinates' parameter. When I use the r.viewshed tool in QGIS and select a coordinate by point and click, the coordinate is (eg.) "3661826.152738,5486199.870317 [EPSG:31467]" and the viewshed analysis works as expected.
However, I'd like to perform a viewshed analysis for a large number of points in a point shapefile, so I used the Add X/Y Field to Layer tool, exported them as csv and am now trying to perform the viewshed analysis for each point using the processing toolbox in Python.
When I manually insert the coordinate "3661826.152738,5486199.870317 [EPSG:31467]", my code works fine. However, the coordinate added with the Add X/Y tool are of the format "11.3445, 49.6174" (numbers here rounded) and the viewshed tool using the coordinate parameter "11.3445, 49.6174 [EPSG:31467]" produce no shapefile (and also no error).
I assume that I have to translate my coordinates from one format into the other and that projections may play a role. How can one do that? I have tried to reproject my raster and project file in [EPSG:31467] as well as in [EPSG:4326], though none of this has changed the coordinate format of add X/Y.


texts:
How to display some information on the label, as an example of the attribute "32710100010080123" too long.

I just want to show four last number to "0123", help me to make vbscript expression label.


texts: How to use ST_Collect to find all polygons with matching ID?

texts:
I have a PostGIS table which holds many polygons. Some polygons have the same ID but are their own discrete area.

How do I do an update query on table which takes groups of polygons with matching ID's and combines them into a single multipolygon geometry.


texts:
I would really like some insight to the pros and cons of one approach versus the other.

I recently wrote some python code that inspects two layers of the same test ESRI GDB file using RTree (bounding box) and Shapely.  This was done on a thin client remotely accessing a Linux server running virtual instance of my Linux distro session (TRANSLATION:  Things run slow).

My code completed evaluating the test GDB (80MB) in ~ 84 seconds.  Very good given the environment I'm working in.

A co-worker of mine did the same evaluation within the same environment by first converting the test GDB file over to a postgreSQL database (~ 16 to 22 seconds) using OGR and then making the same evaluations in about 1 second.

My question to the community given the clear difference in performance is whether or not I should go all in with the PostgreSQL / PostGIS approach or are there cases where Shapely / R-Tree would be better?


texts: Featureclasses with attachments: objects cannot be updated outside edit session

texts:
I am writing an ArcMap add-in using C# with Visual Studio 2010 and ArcGIS 10.1.
I want to add features to a featureclass in an Enterprise Geodatabase (ArcGIS Server 10.1 on SQL Server). This is my code:

var mapLayers = GetAllLayersFromMap();
             if (mapLayers.Count > 0)
             {

                 IWorkspace workspace = ((IDataset)mapLayers[0]).Workspace;
                 IFeatureWorkspace fWorkspace = workspace as IFeatureWorkspace;
                 IWorkspaceEdit workspaceEdit = fWorkspace as IWorkspaceEdit;

                 workspaceEdit.StartEditing(false);
                 workspaceEdit.StartEditOperation();

                 ILayer layerPunt = mapLayers[0] as ILayer;
                 IFeatureLayer featureLayer = layerPunt as IFeatureLayer;
                 IFeatureClass punten = featureLayer.FeatureClass;

                 IFeatureBuffer bufferPunten = punten.CreateFeatureBuffer();

                 IFeatureCursor cursorPunten = punten.Insert(true);

                 workspaceEdit.StopEditOperation();
                 workspaceEdit.StopEditing(true);
             }

When I load a featureclass that has no attachments it works fine. If I create attachments to the same featureclass and try again, it fails. The code crashes on the line where the insert cursor is created: "Objects in this class cannot be edited outside an edit session". When I delete the attachments from the featureclass it works again.

If I start an edit session from ArcMap through the Editor toolbar, I can edit featureclasses with attachments (btw: this proves that the licensing level is ok). So I tried starting the edit session from code through the IEditor interface. It fails on the same line with a COM Exception.

Ideas anyone?


texts:
Using QGIS I have to merge 416 raster files with elevation data, at the beginning I´m using with only 4 of them to try out:

When using "merge" (Raster > Miscellaneous) the result are a file with only black and white which have not included all the raster files i selected.
With "build virtual raster" (Raster > Miscellaneous) I get the merged file but with a different colour scale. I don´t know if it would affect further analysis.
And with the SAGA tool "Mosaic raster layers" the grid file could not be generated.

I have tried modifying some options without success.

How do I maintain the colour scale when using "build virtual raster"?


texts: Raster Values into Gridded Polygons with Same Cell Size using ArcMap

texts:
I have data in the form of a raster with 1 km x 1 km cells. I need to transfer it into a vector polygon grid layer (USNG) which also has 1 km x km cells, but has divergence around UTM zones and does not match up with the raster.
How does one accurately transform raster to polygon data for the same cell size, but different cell locations?
I am using ArcMap.


texts: Problem with setValue() function in ArcGIS 10

texts:
I'm trying to add data to an attribute table in ArcGIS 10 using the following code:

def make_floor_no( shapefile ):
    "Makes header for number of Floors (FLO) and calculates value"

    fieldName = "FLO"
    try:
        ARCPY.AddField_management(shapefile, fieldName, "DOUBLE")
    except:
        print "Field already there"

    # loop through attribute table
    Rows = ARCPY.SearchCursor( shapefile )

    for row in Rows:
        floors = round( row.getValue( 'HGT' ) / 3.0, 0)
        print str(floors)
        row.setValue( fieldName, floors )
        Rows.updateRow(row)

However, I keep getting an error along the line

row.setValue( fieldName, floors )

I can't spot anything wrong with this, and have tried a few different options. Is there something wrong with my syntax?

The error message is as follows:

 File "Z:\ConstructionMaterial.py", line 100, in <module> make_floor_no( InputFile )
 File "Z:\ConstructionMaterial.py", line 71, in make_floor_no
    row.setValue( fieldName, floors )
 File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\arcobjects.py", line 941, in setValue
    return convertArcObjectToPythonObject(self._arc_object.SetValue(*gp_fixargs(args)))
RuntimeError: ERROR 999999: Error executing function.


texts: I want to find base station coverage for a instant location. I have building data. And i determine optimum coverage as 500m far away from a base station (buffer). So how can i find optimum base station location for this example. Looks like that


texts: Building web application to get PostGIS query from user and view result map using GeoServer via OpenLayers?

texts:
I'm a newbie in web-map application area.

I just found the small application in this tutorial interesting

I am a bit confused, because when I publish a map using GeoServer and OpenLayers, I have to create a layer in GeoServer first and this layer is fixed with a table from PostGIS and then use OpenLayers to call GeoServer to render the map.

How can the tutor of that tutorial change the data source of the layer depending on the query input by the users?
It means there are ways to change the data source of one layer on the fly right?

Please advice me how or can anyone tell me where can I download the code of the app in that tutorial.


texts: Reference to Basemap instance's coordinate system

texts:
I use the matplotlib basemap toolkit (mpl_basemap.Basemap) to plot GeoTIFF images on a map background using pcolor() or pcolormesh(). For this purpose, I need to transform a meshgrid of x and y coordinates of the raster data into the coordinate system of the plot, for example polar stereographic. Currently, this is how I accomplish this:

Create an osr.SpatialReference() object and use it to grab the coordinates of the GeoTiff image via its ImportFromWkt() method, called with the GeoTIFF's projection as an argument.
Export this spatial reference to proj4 via .ExportToProj4() and use it to create a pyproj.Proj() object.
Call this class instance with inverse=True on a meshgrid of easting and northing values of the GeoTIFF image to create a lat/long meshgrid
Call the Basemap class instance on this lat/long meshgrid to create a meshgrid in the Basemap object's coordinate system

Is there an easier way that doesn't proceed through a lat/long representation? If I could get a reference to the coordinate system of the Basemap instance, I could set up a coordinate transform directly between the GeoTIFF dataset and the map.


texts: Using stat_density2d with coord_map in ggplot2

texts:
In my example I am plotting a shapefile for all England & Wales, however I'm interested in the greater London area. If I allow ggplot to simply default and display the entire shapefile then the heatmapping generated by stat_density2d isn't very useful for me, I need to change the scale to show the Greater London area.

To address this I need to set the xlim and ylim to effectively zoom in. I can do this in one of two ways:

Use the xlim & ylim standalone functions in ggplot
Use xlim & ylim within coord_map

The former causes the heatmapping to scale correctly (ie size of the tiles reduce to give a meaningful indication of density in a smaller area) but the shapefile data gets truncated to the plot leading to some very odd looking polygons.

The latter causes the heatmapping to remain at the scale of the whole country (ie the tiles are massive and so not useful for showing activity density when the map is zoomed in) but the polygons display correctly.

Is there a way to get ggplot to change the xlim/ylim values of a plot without truncating the shapefile/polygon data so that it also rescales the tiles generated by geom_density2d?

I know I could fudge this by only plotting a subset of the shapefile polygons but this seems...messy. It would be way cooler if I could choose anywhere in the country to zoom/centre on and have everything scale correctly without having to chop and change the shapefile data.

This is a follow up to my previous question: How do I create a heatmap in ggplot?

Edit 1:
I'm currently investigating a few undocumented arguments that stat_density2d passes to the kde2d 2d density function from the MASS package. There are a few of these but I get the impression a combination of n and lims may help. Basically if I pass a higher value to n (this dictates the number of rows and columns there are in the matrix produced) then the 'resolution' of the 2d density estimation increases at the expense of processing time (this is applied to the entire bounding area of the shapefile rather than just the window I'm limiting the plot to using coord_map). I then use lims to try and limit the matrix to just the plotting area but that doesn't seem to have much effect.


texts: How to map the frequency of overlapping points within an urban area?

texts:
I've been asked to do a favor for someone since I'm "the mapping guy" and could use some theoretical help by way of examples since my actual cartography skills are a little rusty.

I have a series of points (~200) that represent events within a city over the course of a year. These events occur all over, but it's common that some points will overlap (e.g. events occurred at the same location).

I'd like to display all of the points, but in the case of overlapping points would like to highlight the fact that the frequency of events is higher at those locations. I can think of a few ways to do this (size, color), but could use some examples or information on techniques that I could use to turn this data into an informative map.

I'm currently using ArcGIS 10 to project and display the points.


texts:
My Mac OX 9.5 freezes or does not respond every time I try to add a CSV file. I am using the following tutorial and following every step, but my computer hates it.

http://fredgibbs.net/tutorials/qgis/linking-and-styling-data-with-qgis.html

as soon as I hit the browse file and select the CSV file that I'd like to import


texts:
I am new to QGIS and trying to understand it to draw historical maps. I run Windows 10 and QGIS 3.10 5. Previously I have used the georeferencer to load map Jpegs but trying to load one yesterday the Georeferencer just didn't open.
I tried it several times, closed the program and windows down and restarted several times. On one occasion and the plug in disappeared from the plugin drop down.I installed the more recent release but the same issue it wont run. I have gone back to 3.10.5 the plugin is there but will not run.


texts: Restrict editing access to my PostgreSQL/ PostGIS database used through QGIS

texts:
We have a small network of users accessing our QGIS data stored on a PostGIS database which sits on a 'Q' drive in my computer. This drive is then shared with the others.
Due to the other users being beginners in the GIS world I created a number of layer files (.qlr) for them which point back to the data in my PostGIS database. They can only see these layer files.

But I've just realised when they use these layer files they can still edit, add and remove the data.

How do I restrict this?

Is there a simple way in QGIS or even Windows? (setting permissions on the sharing folders does not seem to work)

The database was created through something called 'pgAdmin 4' and I've tried adding another user/role with different privileges but I can't seem to get that to work its seems to deny access completely or its all still wide open.

The privileges I can alter seem to be too high up the hierarchy all I need to restrict is the ability to write and delete to the schema's and tables, leaving myself with editing rights.

I have no idea how to write in SQL, I'm a geographer not IT


texts: Why is scale of image different in QGIS 3.10 than created?. CAD import and twf world file seems correct

texts:
I created an image by taking a picture of a tape measure(ruler).  I rotated this image to north and made a world file .tfw for the tiff image.  to get the pixel width for the twf file, I divided the measured width (1.308 meters) by the pixels in the image 1920.  This resulted in 0.00068125 m/pix. I used pyproj with EPSG:3857 to get the upper left corner of the image from the GPS location of where the image was created.  When I bring this into QGIS with projection 3857 it measures smaller by a factor of 0.705.   When I bring it into a cad system, it is correct.
I am using QGIS version 3.10.3-A Coruña on a Windows system. EPSG:3857.
Is there something I need to turn on or scale with QGIS?
EDIT
By just scaling the pixel width number by the cosine of the latitude (from the suggestion below) in the world file this seems to solve the issue for now.


texts:
I'm having issues using OGR to interact with tables on my CartoDB account.  Initially I was trying to upload data using ogr2ogr:

ogr2ogr --config CARTODB_API_KEY abcdefghi -t_srs EPSG:4326 -f CartoDB "CartoDB:thementat" agg.json

This error out, suggesting that I was attempting to create a database... something the CartoDB driver doesn't support:

ERROR 1: HTTP error code : 404
ERROR 1: HTML error page returned by server
ERROR 1: CartoDB driver doesn't support database creation.
ERROR 1: CartoDB driver failed to create CartoDB:thementat

Research had indicated that that I probably shouldn't trust this  so I tried to simplify using ogrinfo:

ogrinfo -ro "CartoDB:thementat tables=parcel"

This returned the following error:

ERROR 1: HTTP error code : 404
ERROR 1: HTML error page returned by server
FAILURE:
Unable to open datasource `CartoDB:thementat tables=parcel' with the following drivers.
  -> PCIDSK
  -> netCDF
  -> PDF
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> OGR_SDTS
  -> S57
  -> DGN
  -> OGR_VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> Interlis 1
  -> Interlis 2
  -> OGR_GMT
  -> GPKG
  -> SQLite
  -> WAsP
  -> PostgreSQL
  -> OpenFileGDB
  -> XPlane
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDUMP
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> OGR_PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> EDIGEO
  -> GFT
  -> SVG
  -> CouchDB
  -> Cloudant
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> CartoDB
  -> AmigoCloud
  -> SXF
  -> Selafin
  -> JML
  -> PLSCENES
  -> CSW
  -> VDV
  -> TIGER
  -> AVCBin
  -> AVCE00
  -> HTTP

I suspect it isn't connecting to my CartoDB account, because this works:

ogrinfo -ro "Carto:gdalautotest2 tables=tm_world_borders_simpl_0_3"

Any ideas?


texts:
I've got some XML documents which conform to a known schema which include geometries in GML format.

I'm looking to perform validation on the XML using XSD and Schematron, but I'll need some way of performing spatial queries within the Xpath language (I presume via extension functions).

I was wondering if anyone is aware of a standard for implementation I can use, or indeed if someone has already done this - I've come up empty on google.

As an example (representative only, only attempting to demonstrate the xpath part of the question (which is the question really - the fact I'm aiming to use it in schematron is moot))

My XML:

<Things>
  <Thing type="A">
    <Geometry>...GML...</Geometry>
  </Thing>
  <Thing type="B">
    <Geometry>...GML...</Geometry>
  </Thing>
</Things>

Xpath to return things of type A which spatially intersect with things of type B (again, I'm making up a function extension namespace and a (pretty dumb) function to give an example of what I'm trying to accomplish):

/Things/Thing[@type='A' and geo:has-intersection(Geometry, /Things/Thing[@type='B']/Geometry)]

As this seems somewhere between development and GIS, I've cross posted on GIS and StackOverflow.


texts:
I try to use a WFS service with OpenStreetMap to request GeoServer GeoJSON file. First I had a error message on the Firefox Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://localhost:8080/geoserver/robot/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=robot:geofield&outputFormat=aoolication%2Fjson. (Reason: CORS header 'Access-Control-Allow-Origin' missing).

Then I pasted the following lines in the web.xml from /usr/share/geoserver/webapps/geoserver/WEB-INF

<filter>
  <filter-name>CorsFilter</filter-name>
  <filter-class>org.apache.catalina.filters.CorsFilter</filter-class>
  <init-param>
    <param-name>cors.allowed.origins</param-name>
    <param-value>*</param-value>
  </init-param>
</filter>
<filter-mapping>
  <filter-name>CorsFilter</filter-name>
  <url-pattern>/*</url-pattern>
</filter-mapping>

After when I go on http://localhost:8080/geoserver/ I have a HTTP ERROR: 503 and still have the same CORS header problem

I don't know what is wrong, maybe this method doesn't work with latest GeoServer versions

If someone has a good and clear answer, it will be great because I think I'm not the only one who want to use a wfs with GeoServer and OpenLayers.


texts: How can I create a button to add a POI layer?

texts: I am making a map of a city with openlayers. I have the shapefile data of the POI (hospital,bank,school etc.) I want to add the layers with buttons not in layerswitcher. Like whereis.com's  quick search buttons ... Can someone please help?


texts: I have two geometries that should be identical but they're not. Where is the difference?


texts: Connecting points with line using virtual layer in QGIS

texts:
I have 10 points in QGIS layer as shown in Image 1.

I am trying to connect points with lines in required sequence as shown in Image 2. However, when I tried to connect all points with virtual layer, the connections are not according to my required sequence (yellow lines) as shown in Image 3. What should I do to get the required connections?

The features in attribute table are shown here. Do I need to change the IDs of points to achieve the required connections?

I was thinking about parent layer and child layer for required connections but I am not sure how to achieve that.

The shapfile is available here

Virtual layer query that I have tried (Image 3):
In this query, I am trying to connect points with line by saving the start and end points of each connecting line.
with
build_segment as (
select
  a.GlobalID as start_point,
  b.GlobalID as end_point,
  makeline(a.geometry, b.geometry) as geometry
from Points_Colac_Demo_advanced a
left join Points_Colac_Demo_advanced b on (a.OBJECTID + 1 = b.OBJECTID))

select
  *
from build_segment
where geometry is not null


texts: Geoserver - generating tiles - bbox is incorrect

texts:
I am beginner in geoserver. I would like to generate tiles in geoserver. I have added the layer to geoserver. When i choose the option to 'Preview layer' i can see the layer with the correct bbox:

Then, i go to gwc, i choose the option 'seed this layer', and feel in this part as follows:
Number of tasks to use: 01
Type of operation:  Seed - Generate missing tiles
Grid Set:   EPSG:4326
Format: image / png
Zoom start: 00
Zoom stop:  02
Bounding box:   -4.601615515076983 39.8769407866263 -3.0527184873294764 41.16710732525929

And as a result i can see almost all of the tiles completely white and one tile with black point...

Why does it happens?


texts: Exporting Feature Class to Shapefile in ModelBuilder?

texts:
ArcGis 10.4

I have to create buffers in the ModelBuilder.
At one point I have a feature class from which I have to spatially select some items (via tool "select by location"). But this tool does not allow a feature class as an input parameter.

The only thing that works is to seperate the work flow. If I run the model until this point and then manually take the produced output, the select by location tool accepts this as an input. But of course that destroys the workflow and will not let me do any iterations for example.

Unfortunately I cant find an option to export this featureclass (to a shapefile?) and google is not helping me either.
Any ideas?


texts:
i´m trying to load the layer Mroad_Lden_17 into R with the code below:
proj_folder <- "C:/Users/... "
setwd(proj_folder)

# helper function for querying uba
load_laerm <- function(what) {
    filename <- file.path(proj_folder, paste0(what, ".RData"))
    if(!file.exists(filename)) {
        message("building ", filename)
        wfs_laerm <- "https://datahub.uba.de/server/services/INSPIRE/end/MapServer/WFSServer?"
        laerm_client <- ows4R::WFSClient$new(wfs_laerm, serviceVersion = "2.0.0")
        laerm_url <- httr::parse_url(laerm_client$url)
        laerm_url$query <- list(
            service = "wfs", version = "2.0.0", request = "GetFeature",
            typename = paste0("INSPIRE_end:", what), srsName = "EPSG:3857"
        )
        request <- httr::build_url(laerm_url)
        ans <- sf::read_sf(request)
        save(ans, file=filename)
    } else {
        message("loading from ", basename(filename))
        load(filename)
    }
    stopifnot(inherits(ans, "sf"),  inherits(ans, "data.frame"))
    message("==> ", what, " loaded.")
    return(ans)
}
Mroad_Lden_17 <- load_laerm("Mroad_Lden_17")

Unfortunately, the layer is really big, and R needs much time to load it, e.g. it stops the loading process, because the layer is too big with its 809559 features. I thought about a way to load the layer in "smaller pieces" and put it togther afterwards, but I do not know how to do that.


texts: org.geotools.jdbc.JDBCDataStore: add information that is stored with the shape file?

texts: Is there a way to add metadata such as file name and file creation timestamp as separate table columns in the table when org.geotools.jdbc.JDBCDataStore stores a shape file in a PostgreSQL database?


texts: Delete all points outside the contour of a WMS raster map

texts:
I have two layers in QGIS. One WMS raster layer with a country map and one layer with points. It looks like this:

I want to remove all the points that lie outside the contours of the country map, so I can just keep the points that lies inside the country. How can I do this easier than removing all the points one by one?


texts: Why do I lose data after converting from .xlsx to .shp?

texts: I have a fairly large excel file that I need to convert to a shapefile, but everytime I export the data, I lose a bunch of it. The excel file has about 277,000 records, but the shapefile will only contain about 260,000 records. I've tried to convert the .xlsx to .csv, which only lost about 2,000 records the first time I exported it, but lost about 20,000 the second time. Why am I losing data during this conversion?


texts: Automatically apply topology rule to only selected portions of a .shp file?

texts:
I am looking for some examples/tutorials on how to apply a topology rule to only a selected portion of a shp file using a script in ArcGIS/GRASS.

Using ArcGIS seems to always require an interactive approach - with the topology editor able to highlight exeption and the user can mark them as exeptions. GRASS seems to have much less detailed topology editing capabilities compared to ArcGIS.

I would like to have a script that automatically fixes exceptions or not depending on something like an sql statement or the proximity of exceptions to another vector file. For example, I want to fix dangles in roads everywhere but downtown, because downtown they are purposeful but out of downtown they are mistakes.


texts: Create two new polygons inside existing polygon hole

texts:
I would like to create new polygon inside the existing polygon.
It does not means donuts polygon.

I would like to have two polygon without overlapping.

I can create if I don't check "Avoid intersection" from snapping option.
Unfortunately, there is overlapping new polygon and existing polygon.

It must not have overlapping area.


texts: Rendering thumbnails of a GPX track with my own tiles as base layer

texts:
I created my own Map using Tilemill. Now I'd like to display thumbnails of a GPX track where my map is the base layer.

Do you have any idea how to do this performantly? I wrote a script that does it using Tilemill but doing that for a couple thousand tracks would take forever...

I'd have all tiles pre-rendered anyway as either mbtiles or just as tiled pngs in folders...

Any ideas are very welcome :)

:: EDIT ::

If you could tell me a tool that takes a datasource (mbtiles or web) like https://github.com/vvoovv/tile-utils for example but that would be able to clip the result to the bounding box i specify rather than taking all the tiles that fit into the bounding box, I should be all set


texts:
I have a script that runs through a database and deletes all the empty feature classes.  Often I deal with databases that have hundreds of empty feature classes and it can take quite a while to run the tool that deletes the empty FCs.  I have an idea to delete the entire feature data set that houses the empty FCs instead (thinking it would be faster).  Any suggestions on how to do this?

Essentially, if a feature data set has a FC in it that has data it would still have to individually delete the empty FCs, but if the entire feature data set only houses empty FCs, I would like to just delete the data set itself.  It may actually not be faster (just an idea). Here is my script as it stands.

import arcpy
arcpy.env.workspace = arcpy.GetParameterAsText(0)

datasetList = arcpy.ListDatasets('*', 'Feature')
for dataset in datasetList:
    fcList = arcpy.ListFeatureClasses('', '', dataset)
    for fc in fcList:
        results = arcpy.GetCount_management(fc)
        count = int(results.getOutput(0))
        if count == 0:
            arcpy.Delete_management(fc)
            arcpy.AddMessage('Deleting ' + fc)
        else:
            pass