texts: Why is r.out.gdal exporting empty rasters?

texts:
I need to export a raster I produced with GRASS. The Map Display portraits this raster like this:

The metadata of this raster also looks pretty normal:

> r.info map=blenddem
 +----------------------------------------------------------------------------+
 | Map:      blenddem                       Date: Fri May  5 11:35:05 2017    |
 | Mapset:   DEMs                           Login of Creator: desouslu        |
 | Location: Luzern                                                           |
 | DataBase: /home/desouslu/maps/grass              |
 | Title:                                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    DCELL                                                      |
 |   Rows:         3369                                                       |
 |   Columns:      3110                                                       |
 |   Total Cells:  10477590                                                   |
 |        Projection: CH1903+ / LV95                                          |
 |            N:   210553.5    S:     208869   Res:   0.5                     |
 |            E:     668202    W:     666647   Res:   0.5                     |
 |   Range of data:    min = 430.691497802734  max = 600.981266929298         |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.patch                                                    |
 |                                                                            |
 |   Comments:                                                                |
 |    r.patch input="tmp_3151920,uavdem" output="blenddem"                    |
 |                                                                            |
 +----------------------------------------------------------------------------+

I am using r.out.gdal to export this raster, and apparently everything goes on fine:

> r.out.gdal input=blenddem output=/home/desouslu/maps/blenddem.asc format=AAIGrid --overwrite -c
Driver <AAIGrid> does not support direct writing. Using MEM driver for
intermediate dataset.
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float64>
Input raster map contains cells with NULL-value (no-data). The value -nan
will be used to represent no-data values in the input map. You can specify
a nodata value with the nodata option.
Exporting raster data to AAIGrid format...
 100%
r.out.gdal complete. File
</home/desouslu/maps/blenddem.asc> created.

However, the resulting raster is actually empty:

> cat /home/desouslu/maps/blenddem.asc
ncols        1
nrows        1
xllcorner    0.000000000000
yllcorner    0.000000000000
cellsize     1.000000000000
NODATA_value  -nan
 -nan

What is exactly going on here? Am I using r.out.gdal incorrectly?


texts: Warp polygons/lines from distortion grid, or similar

texts:
Anyone have any ideas how I might go about warping polygons or line based on point based distortion map?

Specifically, let's say I have a (non-regular) point data set. I create a clone of the data set and translate points in a non-regular fashion. Is there a way to use the difference of the two layers to warp a related line or polygon dataset? It doesn't need to be "accurate" (whatever that means in this instance), just give a rough indication of what the underlying topology looks like giving the applied distortion.

Using either PostGIS, QGIS, Shapely, GEOS etc. Or if you have an Arc* version I can't probably make that work too.


texts:
My specifications :

Python 3.9.1
gdal 3.1.4
rasterio 1.1.8
Windows 10 Business
conda version 4.9.2

I have a raster that has no crs :
os.chdir(path)
file_list= os.listdir()
file_list[0]
>>> RGEALTI_FXX_0605_6670_MNT_LAMB93_IGN69.asc

with rasterio.open(liste_fichiers[0]) as src:
    print(src.profile)
    print(src.crs)
    print(src.bounds)
    print(src.indexes)
>>> {'driver': 'AAIGrid', 'dtype': 'float32', 'nodata': -99999.0, 'width': 1000, 'height': 1000, 'count': 1, 'crs': None, 'transform': Affine(5.0, 0.0, 604997.5,
       0.0, -5.0, 6670002.5), 'tiled': False}
None
BoundingBox(left=604997.5, bottom=6665002.5, right=609997.5, top=6670002.5)
(1,)

For your information, I do have a GDAL_DATA environment variable. (I found this information here).
When I open the raster with QGIS, I can in fact notice that there is no CRS, but on the other hand I can easily select one.
I am looking forward to do the same with rasterio. Here is the code I tried to do so :
with rasterio.open(liste_fichiers[0], mode='r+') as src:
    print(src.profile)
    src.crs = rasterio.crs.CRS.from_epsg(2154)
>>> {'driver': 'AAIGrid', 'dtype': 'float32', 'nodata': -99999.0, 'width': 1000, 'height': 1000, 'count': 1, 'crs': None, 'transform': Affine(5.0, 0.0, 604997.5,
       0.0, -5.0, 6670002.5), 'tiled': False}

---------------------------------------------------------------------------
CPLE_AppDefinedError                      Traceback (most recent call last)
<ipython-input-191-a0546e13d796> in <module>
      1 with rasterio.open(liste_fichiers[0], mode='r+') as src:
      2     print(src.profile)
----> 3     src.crs = rasterio.crs.CRS.from_epsg(2154)

rasterio\_base.pyx in rasterio._base.DatasetBase.__exit__()

rasterio\_base.pyx in rasterio._base.DatasetBase.close()

rasterio\_io.pyx in rasterio._io.BufferedDatasetWriterBase.stop()

rasterio\_io.pyx in rasterio._io._delete_dataset_if_exists()

rasterio\_err.pyx in rasterio._err.exc_wrap_int()

CPLE_AppDefinedError: Deleting RGEALTI_FXX_0605_6670_MNT_LAMB93_IGN69.asc failed: Permission denied


texts: updateCursor not updating rows correctly with conditional statement only with arcpy toolbox

texts:
I am attempting to create a script with arc toolbox in which part of the script populates a probability field based on thresholds on a percentage field.  I would like the thresholds for the ranges upon which to classify the probability field to be parameters that the user can input in the tool.  However, when I attempt to run the tool, it does not properly update the fields.  Instead it fills with the last value (in the case of my script below, it populates everything with "L").  The script runs fine (through toolbox or stand-alone) if I replace the parameters in the script with values. I have tried defining these threshold parameters as "double" or "any value" in arc toolbox, both with no luck.  Does anyone know why I am having this problem?

# populate occurrence probability based on percentage thresholds
with arcpy.da.UpdateCursor(zoneTable, ["PERCENTAGE", "probability"]) as cursor:
    for row in cursor:
        if (row[0] >= highThreshold):
            row[1] = "H"
        elif (row[0] < highThreshold and row[0] >= medThreshold):
            row[1] = "M"
        else:
            row[1] = "L"
        cursor.updateRow(row)

highThreshold = arcpy.GetParameterAsText(2)
medThreshold = arcpy.GetParameterAsText(3)


texts: Is there a minimap plugin(function) in QGIS?

texts:
I'd like to add minimap(or mini-map) feature to my app, which were developed based on QGIS. I know there is a overview feature, however, it is quite different from minimap. The minimap should show basic features of canvas window, such as map, item.
It should be like this:


texts: Representation of spatial lines with altitude

texts:
Is there any R data structure that supports spatial lines that include height information? I'm looking for something like SpatialLinesDataFrame with support for altitude in the lines. Associating a time with each point would be nice but not essential.

My data is airplane flight tracks. I want to be able to do spatial analysis such as finding the ground distance and altitude at the closest point of approach to a fixed point on the ground, and visualizations such as maps where the segments are coloured by altitude.

The data is supplied as KML files. Currently I am using rgdal::readOGR to read the lat/lon data and metadata, and maptools::getKMLcoordinates to read the same data as a data.frame containing lat/lon/altitude. This works but it is awkward to use two data structures to contain pieces of the same data.

Any suggestions?

(I originally posted this on StackOverflow but no answers there...)


texts:
I have two Landsat 5 NDVI images (Surface Reflectance products) from the same period of the year (only a few days apart), but one is from the 1980s and one from the 2000s. I am trying to do a change detection procedure based on NDVI differencing, and the expected outcome (based on aerial photographs) is that NDVI should have increased due to a widespread reforestation process, so the result of the differencing process (2007NDVI - 1987NDVI) should highlight mainly a) unchanged areas and b) areas of ndvi increase.

I tried to reclassify the difference image using a threshold calculated as μ ± n·σ; (μ represents the mean, and σ the standard deviation). I used n=1.
(reference paper: http://www.sisef.it/iforest/contents/?id=ifor0909-007 )
The difference image classified this way, however, shows mainly unchanged pixels, so is it possible that the difference image shows no increase simply because the data are not "normalized"? Since this is a SR product, it has already been atmospherically corrected, so it may be that due to climatic differences during those two years in one case the vegetative season was more advanced than the other.

These are the image histograms:

The first is from 1989 and the second from 2007.

Would in this case a histogram matching procedure be of any help (or any other procedure you might suggest)? And if so, how can I perform that using QGIS? I couldn't find any dedicated tool or plugin for it.


texts: How to get LatLng array from geoJSON Layer in Leaflet?

texts:
i need to hand the latlng coordinates from my phpPostGIS-served geoJSONLayer as an array over to a new added HeatMap Layer. L.geoJSON has an build-in option: coordsToLatLng(  coords )" and a static method for that: coordsToLatLngs(  coords,  levelsDeep?,  reverse? ).
I don´t get any errors, but the heat map does´t show up. If i build a new L.heatLayer(coord).addTo(map); outside the function, it works…
Does anyone have an idea?
This is what i have so far:

    $.getJSON("postgis_geojson.php?geotable=shops&geomfield=the_geom", function(data) {

        L.geoJson(data, {
            pointToLayer: function (feature, latlng) {
                return L.circleMarker(latlng, geojsonMarkerOptions);

            },
             onEachFeature: function (feature, layer) {
                 popupOptions = {maxWidth: 200};
                layer.bindPopup(feature.properties.popupContent);
            },
            coordsToLatLng: function (coords) {
            return new L.heatLayer(coords);
            }

        }).addTo(map);
        });


texts:
I'm trying to merge 9 raster files I have. for some reason. if I merge a few of them it works, especially when they are close to each other. But whenever I try to merge all of them, together I get this bad result:

I assume that maybe it has to do with the fact that some of the rasters are very far from each other geographicaly but I don't know if it has to do with that.

My end goal is to merge all the 9 rasters into one raster layer.


texts: Aligning environmental layers causing a change in nodata values

texts:
I'm trying to align and re-sample environmental layers in QGIS (2.18.3). The layers I am using were obtained from the world clim data base Bio3,Bio18,Bio12,Bio10 as well as a Sheep and Cattle layer obtained from the FAO website. The World clim layers are at pixel size of 0.0416667,-0.0416667 and the two livestock layers are at a pixel size of 0.05,-0.05. Using the identify tool in QGIS I can see the ocean areas in this raster have a value of 0 for the livestock layers and a nodata value in the environmental layers (note my analysis is only interested on the land). I have been trying to use the raster align tool in QGIS so I can match these layers for some later Maxent analysis. I have used the following procedure in QGIS. Raster>Align Raster. Then I have been adding all the 6 layers, making the livestock layer the reference layers as it is of lower resolution. Checking the CRS value is at WGS 85, as well as adding coordinates to trim. This process occurs correctly, with pixel size matched. However my concern arises when I look at the trim and aligned rasters. Three raster Bio 10,12,18 all get a minimum value of -32768, were as Bio3 gets 43. When I used the identify tool again on the rasters Bio10,12,18 in the oceans show a value of -32768, whereas Bio3 give a value of 255 and the livestock layers are at 0.

Should the min and max values be matching in these layers if I am to compare between them?

Or are these values just proportional?


texts: It is possible to use a GDAL in plugin qgis?

texts:
I work on qgis plugin, and I say why not not work with GDAL, it'll be too useful.
it is possible to integrate the GDAL library to develop QGIS plugin?


texts:
I am trying to zoom to a point in ArcGIS JavaScript API (3.9). I am doing this one of two ways, I am either grabbing a point feature from ArcGIS Server (10.1), or I am grabbing the center of a polygon and zooming to that.

I was getting an error "exception in animation handler for: onAnimate" when zooming to any point. I found here that adding a setTimeout() around your map.centerAndZoom() solves the problem. However, that was only for the point feature being brought back from the server.

The problem still persists when I try to get the center point of a polygon:

I have tried both increasing and decreasing the timeout time, but I have not had any success. The map will seem to zoom to the location, but then immediate zoom back to its original location and then spit this error:

I have found a little research done on "Cannot read property 'left' of undefined" and the most I could find is something about symbology - So, just to help out, here is the symbology as it is in ArcGIS Server:

Renderer:
  Simple Renderer:
  Symbol:
           Style: esriSFSSolid
           Color: [0, 0, 0, 0]
           Outline:
                   Style: esriSLSSolid
                   Color: [76, 230, 0, 255]
                   Width: 1
 Label:
 Description:
 Transparency: 0
 Labeling Info:

If anyone could shed some light on this problem, I would be grateful.


texts:
I have a number of KML files and I would like to convert them to Shapefile format.

Is it possible to build a model that iterates over a number of KML files?

I tried to use File iterator, but it doesn't seem to work.


texts: How to use the gdal.Warp function with options?

texts:
I try to use the gdal.Warp function with options :
# align GFC tree cover with segments
glad_raw = pm.getDwnDir() + aoi_name + '_treecover2010.tif'
glad_aligned = pm.getDwnDir() + aoi_name + '_treecover2010_aligned.tif'

if not os.path.isfile(mask):
    return 'No treecover file'

options = gdal.WarpOptions(
    xRes             = res['x'],
    yRes             = res['y'],
    creationOptions  = "COMPRESS=LZW",
    outputBounds     = [int(v) for v,i in enumerate(bb)],
    dstSRS           = proj
    #options          = 'overwrite'
)
ds = gdal.Warp(glad_aligned, glad_raw, options)
ds = None

but I get the following error :

TypeError: Warp() takes 2 positional arguments but 3 were given

strange because the help(gdal.warp) gave me Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)


texts: QGIS wrong measurments on Google satellite images

texts:
I am new to QGIS but I understand the basic principles. I'm working with QGIS 3.8.3 on OpenSUSE Tumbleweed Linux using the KDE environment.

Still trying to figure out QGIS, I'm opening a brand new project with default settings. I load Google satellite images in a layer and start measuring an area that I own and have exact measurements of. A 200 meter tree line is reported by QGIS to be 230 m, and another tree line 90 degrees to the first one is 210 m is reported by QGIS to be 240 m. These measurements are correct when using Google maps, so it's a QGIS setting issue.

Searching for answers I made sure that the project projection is WGS pseudo Mercator 3857 as this is supposed to be the projection of Google maps but it made no difference.

Any suggestions on how to set this up correctly?


texts: Stratified sampling in Google Earth Engine

texts:
I am trying to get stratified random points out of GHSL (Global Human Settlement Layer) which is a classified raster over my study area (partially Europe and Asia continent), see the picture below.

I face multiple errors like computation time out.

Here is the code:

var ghsl  = ee.Image('JRC/GHSL/P2016/BUILT_LDSMT_GLOBE_V1')
var ghsl1 = ghsl.clip(aoi).select('built')

var masked = ghsl1.updateMask((ghsl1.eq(1)).or(ghsl1.eq(2)).or(ghsl1.eq(5)).or(ghsl1.eq(6)));
var ghsl2 = masked.remap([1,2,5,6],[0,1,2,2])
print(masked)

var ghsl5 = ee.Image(masked);
print(ghsl5)

print('Projection, crs, and crs_transform:', ghsl5.projection());

var stratified = ghsl5.addBands(ee.Image.pixelLonLat())
    .stratifiedSample({
      numPoints: 1000,
      classBand: 'built',
      projection: 'EPSG:3857',
      scale: 38,
      region: features.geometry()
    }).map(function(f) {
      return f.setGeometry(ee.Geometry.Point([f.get('longitude'), f.get('latitude')]))
    })

print (stratified.reduceColumns(ee.Reducer.frequencyHistogram(),['class']));
Map.addLayer(stratified);


texts: Finding points within elliptical area using PostGIS?

texts:
Using PostGIS 2.0.0 I need to filter the results of a function, which returns many POINTs, to find the ones that lie within an elliptical area (like a great circle line with a buffer around it) defined by two points and a distance (all parameters passed into the query).

A filter using ST_Distance does the trick, but I wonder if there's any way to speed it up?

Simplified query:

SELECT candidate_point
FROM my_function(some, params)
WHERE
    ST_Distance(candidate_point, 'SRID=4326;POINT(1 2)'::geography) +
    ST_Distance(candidate_point, 'SRID=4326;POINT(3 4)'::geography) <= 60000

I tried using ST_DWithin as an additional filter, but that slows it down tenfold! Because candidate_point comes from a function, not a table, there is no index to use.


texts: Test Schema Lock does not catch open feature class from within ArcGIS Pro

texts:
I am working on a script tool to be run in ArcGIS Pro (version 2.4.0 with Python 3.6.8) and wish to give the user a warning message and also time to close a feature class so as to avoid the script breaking with a schema lock run time error. I tried applying the code from Pause ArcPy script until feature class lock goes away, however I find that the open feature class is caught when the code is run from an IDE (I'm using Spyder), but it is not caught when run from ArcGIS Pro's Interactive Python Interpreter window or when included in a script tool.

I've tried to make a minimum reproducible example (below). After making the test dataset, Example_Polygons, I open it in Pro and run the test schema lock code from both an IDE and the Pro Python window. lockTest is False in the IDE and True in Pro.

import arcpy, time

# Create test polygon for minimal reproducible example

arcpy.env.overwriteOutput = True
if not arcpy.Exists(r"C:\Temp\test.gdb"):
    arcpy.management.CreateFileGDB(r"C:\Temp","test.gdb")
arcpy.env.workspace = r"C:\Temp\Test.gdb"
srWGS84 = arcpy.SpatialReference("WGS 1984")
arcpy.management.CreateFishnet("Example_Polygons",
                               "0 0", "0 1", 1, 1, 4, 6, None,
                               "NO_LABELS", "DEFAULT", "POLYGON")
arcpy.management.DefineProjection("Example_Polygons", srWGS84)

# Code to compare between when run from IDE and from python interactive window (or potentially script tool)

# Test schema lock
lockTest = arcpy.TestSchemaLock('Example_Polygons')
print(lockTest)

# Give user message and time to close feature class
while lockTest != True:
    time.sleep(10)
    lockTest = arcpy.TestSchemaLock('Example_Polygons')
    arcpy.AddMessage('Please close Example_Polygons')

if lockTest == True:
    arcpy.AddField_management(in_table = 'Example_Polygons', field_name = 'wave_count', field_type = 'SHORT')

I came across https://support.esri.com/en/technical-article/000011739, which states that "The 64-bit background processing cannot gather the schema lock information across applications running as 32-bit." However, I don't think that is relevant in my case because I'm not running anything as 32-bit (to my knowledge).


texts: How to refresh the attribute table of a layer after editing the table in the DB-Manager?

texts: Is there any way to refresh the attribute table of a layer after editing the table (e.g. add a column) in the DB-Manager?


texts: Seeking Android library to run Dijkstra on shapefile?

texts:
I have a shapefile and I want to find the best way of using Dijkstra on Android platform.

Is there any library that I can use it?

I have to convert it to graphs and then run algorithm. I am confused how make things together for this purpose. If it is possible please suggest some libraries that I can use it.


texts: How to use GDAL to burn a selected value from one raster into another?

texts:
I think this should be a relatively simple process, but I'm not quite sure how to approach it.

I have a shapefile containing 1000's of polygons. I also have 19 rasters. All of the polygons fall withing the extent of the rasters. The rasters are equal in extent and resolution. Previously I have used GDAL to extract data for each of these polygons from each of the 19 rasters. This worked very well. However, I now wish to extract the pixels into xyz format or similar for all the areas not covered by these polygons.

Currently I am using gdal_rasterize to burn a chosen nodata value into each of the rasters indicating where the polygons occur. I would then need to translate the raster removing all of the nodata values. This is a very heavy process and is taking an awefully long time for 1 raster let alone 19.

My question is - is there a way to gdal_rasterize all of my polygons onto just 1 of the rasters, then use the nodata values from this raster as a mask, and apply these nodata values to the same locations on the other 18 rasters. I'm sure this would be quicker than my current method, but I don't know how to do this in GDAL.

Any help will be very much appreciated!


texts: What is the best way to take a WKB value in SQL Server Management Studio and create a PostGIS geometry in PGAdmin III?

texts: Normally, I would just use WKT strings but I'm investigating some possible geometry weirdness with the Text methods in SQL Server. I would like to copy/paste the WKB values from the SQL Server Management Studio data grid and paste into PGAdminIII  SQL Editor.


texts: Removing polygon sliver other than by using graphic-to-feature

texts:
I often have to edit the edges of polygons to match or be cut by the edges of other polygons with complex shapes. This sometimes leaves disconnected slivers.

Is there a simpler/faster way to remove these slivers (for example, like the "Eraser" tool in Paint)?

Currently, I make a graphic that contains the part I would like to erase, convert it to feature, and use "Clip" from the Edit menu to discard the area that intersects. Then I have to delete the feature layer that I just created for the sole purpose of erasing the sliver. This works, but is a lot of steps when you have to do it often.


texts: Finding values in value list by typing part of the value in QField

texts:
In QField you have to select a values from drop down value list by scrolling to the desired value. When the value list is very long it would be very helpful if you could enter only a part of the value or only the first letter to select the desired value much faster.
In QGIS you can use the keyboard. Is there a similar option in QField available?


texts:
hope everyone's doing well.

I have a statewide point data set (parent table: 600K rows, 500 columns) I need to analyze at multiple scales based on several different spatial divisions. I'll also need breakdowns by a few key attributes. Here's a list of planned child table types with counts:

mainRegions (4 child tables, a few extra columns for demographics)
subRegions (11 child tables, a few extra columns for demographics)
counties (67 child tables, a few extra columns for demographics)
KeyAttributeOne (15 child tables, no additional columns)
KeyAttributeTwo (48 child tables, no additional columns)

Concerning spatial divisions, we aggregate and report exclusively on the region, sub-region and county level of the state (see image) but analyze the data at the census tract or block level. There are only two attributes in the parent table that can be used to group the records in any meaningful way.

It's my understanding that implementing an inheritance strategy will allow for easy aggregation by geometry, attribute, or a combination of both. I wonder, however, if the easy ad-hoc queries will be worth the work required up-front.

PostGIS In Action (awesome resource) gives a good example of a linestring table sorted by state_fips but I haven't heard any firsthand accounts of using inheritance with geometries.

-Rob


texts: How does MOBAC handle groups in customWmsMapSource?

texts:
I'm hoping to get a custom WMS map source added to Mobile Atlas Creator (v1.9.16). This program is normally available at: http://mobac.sourceforge.net/

I plan to use Maine's Statewide Ortho WMS Services URL:
http://mapserver.maine.gov/wms/mapserv.exe?map=c:/wms/orthos.map&

From there, many groups are available, and within each are sundry sub-groups. How does MOBAC handle culling out single sources from groups and sub-groups?

Here's a sketch of what I've started:

<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<customWmsMapSource>
    <name>MeGIS WMS Low_Tide_2008</name>
    <minZoom>0</minZoom>
    <maxZoom>18</maxZoom>
    <tileType>PNG</tileType>
    <version>1.1.1</version>
    <layers>lowtide2008</layers>
    <url>http://mapserver.maine.gov/wms/mapserv.exe?map=c:/wms/orthos.map&amp;</url>
    <coordinatesystem>EPSG:4326</coordinatesystem>
    <aditionalparameters>&amp;EXCEPTIONS=BLANK&amp;Styles=</aditionalparameters>
    <backgroundColor>#000000</backgroundColor>
</customWmsMapSource>

Here's the impotent outcome:


texts: Converting GeoJSON to Shapefile using ogr2ogr fails

texts:
I am trying to convert GeoJSON to Shapefile with an extent. The Shapefile is generated successfully, however it has no extent and data. Without extent and data, SAGA fails to run using the output Shapefile. Below is the ogr2ogr code.
ogr2ogr -f "ESRI Shapefile" -spat -88.987035 25.585476 -88.987035 30.160476 -clipsrc spat_extent "D:/MSU/watermonitor/monitor/temp/GKS7wuW1Rs9FIv-jAAAA/1598937441196.shp" "D:/MSU/watermonitor/monitor/temp/GKS7wuW1Rs9FIv-jAAAA/1598937441196.json"

This is how the source GeoJSON looks like:
{
  "type": "FeatureCollection",
  "bbox": [
    -88.987035,
    25.585476,
    -88.987035,
    30.160476
  ],
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:EPSG::4326"
    }
  },
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -88.987035,
          25.585476
        ]
      },
      "properties": {
        "id": 239,
        "chlorophyll_a": 52,
        "capture_datetime": "2020-08-30T01:01:18.657+00:00",
        "equal": 1,
        "quantile": 1,
        "quantile_range": "51 - 53",
        "equal_range": "51 - 54"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -88.987035,
          28.964476
        ]
      },
      "properties": {
        "id": 346,
        "chlorophyll_a": 53,
        "capture_datetime": "2020-08-30T03:23:48.477+00:00",
        "equal": 1,
        "quantile": 1,
        "quantile_range": "51 - 53",
        "equal_range": "51 - 54"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -88.987035,
          28.964476
        ]
      },
      "properties": {
        "id": 511,
        "chlorophyll_a": 55,
        "capture_datetime": "2020-08-29T17:22:44.598+00:00",
        "equal": 2,
        "quantile": 2,
        "quantile_range": "54 - 57",
        "equal_range": "55 - 57"
      }
    }
  ]
}

The ogrinfo of the file is the following:
        D:\MSU\watermonitor\monitor\temp\GKS7wuW1Rs9FIv-jAAAA>ogrinfo -so -al 1598937441196.shp
    INFO: Open of `1598937441196.shp'
          using driver `ESRI Shapefile' successful.

    Layer name: 1598937441196
    Metadata:
      DBF_DATE_LAST_UPDATE=2020-09-01
    Geometry: Point
    Feature Count: 0
    Extent: (0.000000, 0.000000) - (0.000000, 0.000000)
    Layer SRS WKT:
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]]
    id: Integer (9.0)
    chlorophyl: Integer (9.0)
    capture_da: Date (10.0)
    equal: Integer (9.0)
    quantile: Integer (9.0)
    quantile_r: String (80.0)
    equal_rang: String (80.0)

And this is what I got from QGIS (Information from Provider):

As you can see, there is no extent/bounding box specified. And it has no data it it even though the source data has 11 rows.
Is there any parameter to that fixes this issue?


texts: How to "save" svg symbols when transforming a map in QGIS using qgis2web to create a web map?

texts:
I'm very new to QGIS so might be getting the terminology wrong.

I created a map with several layers in QGIS. I needed different markers than the ones provided so I used svg markers. Now I'd like to use qgis2web to create a web map - when I do so, my svg markers don't show up but only the Labels which I've also put.

How can I fix this so that the svg markers show up in the web map?


texts: Import Error: No module named shapely.wkb

texts:
OpenERP GeoEngine
I have installed OpenERP 6.1 version, works fine. When I'm trying to add CamptoCamp GeoEngine Modules with OpenERP getting below find error:
Installation Prerequisites:
OpenERP 6.1,CamptoCamp GeoEngine modules(sales,partner) ,Python 2.6,Shapley 1.2.18,Geojson 1.0.1,   POSTGIS 1.5

When starting server in service it throws below error:
12:43:00,283 1904 INFO ? openerp: OpenERP version 6.1-20130918-233418

12:43:00,283 1904 INFO ? openerp: addons paths: E:\OpenERP 6.1-20130918-233418\Server\server\openerp\addons

12:43:00,283 1904 INFO ? openerp: database hostname: localhost

12:43:00,283 1904 INFO ? openerp: database port: 5432

12:43:00,283 1904 INFO ? openerp: database user: openpg

12:43:00,315 1904 INFO ? openerp.service.netrpc_server: starting NET-RPC service on 0.0.0.0:8070

12:43:00,315 1904 INFO ? openerp.netsvc: Starting 1 services

12:43:00,361 1904 INFO ? openerp.wsgi.core: HTTP service (werkzeug) running on 0.0.0.0:8069 12:43:00,440 1904 INFO ? openerp.addons.web: embedded mode

12:43:00,642 1904 CRITICAL ? **openerp.modules.module: Couldn't load module web**

12:43:00,642 1904 CRITICAL ? **openerp.modules.module: No module named shapely.wkb**

12:43:00,642 1904 ERROR ? openerp: Failed to load server-wide module `web`.
The `web` module is provided by the addons found in the `openerp-web` project.
Maybe you forgot to add those addons in your addons_path configuration.
Traceback (most recent call last):

  File "openerp-server", line 252, in <module>

  File "E:\OpenERP 6.1\Server\server\.\openerp\modules\module.py", line 419, in load_openerp_module

  File "E:\OpenERP 6.1\Server\server\openerp\addons\web\__init__.py", line 29, in wsgi_postload

  File "E:\OpenERP 6.1\Server\server\openerp\addons\web\common\http.py", line 453, in __init__

  File "E:\OpenERP 6.1\Server\server\openerp\addons\web\common\http.py", line 523, in _load_addons

  File "E:\OpenERP 6.1\Server\server\.\openerp\modules\module.py", line 139, in load_module

  File "E:\OpenERP 6.1\Server\server\openerp\addons\base_geoengine\__init__.py", line 8, in <module>

  File "E:\OpenERP 6.1Server\server\openerp\addons\base_geoengine\geo_model.py", line 11, in <module>

  File "E:\OpenERP 6.1\Server\server\openerp\addons\base_geoengine\geo_field.py", line 9, in <module>

**ImportError: No module named shapely.wkb**

2013-09-24 12:43:00,642 1904 INFO ? openerp: OpenERP server is running, waiting for connections...


texts:
I am copying a tons of features from a file geodatabase to a database server. First I created an enterprise geodatabase connenction and after a database connection. In the first step I added a database name. After I copy all the features from dekstop to this SQL Server Enterprise Geodatabase the features name modified with a prefix, which is the database name.

Someone know why and what is the best solution to delete all of this prefiex?


texts: How to represent bimodal attribute value distributions on a map

texts:
I'm looking for a visualisation method that enables me to represent bimodal or even multimodal attribute value distributions on line features.

A bimodal distribution has two different modes that appear as distinct peaks in the probability density function (image #1). Mean, median and standard deviation values are of little use to describe such data.

#1 Example of a bimodal distribution:

(source: wikimedia.org)

My use case is the following: For all road links in a network, I have associated speed measurements (many hundreds per link) in bimodal - sometimes multimodal (three modes) - distributions. Now, I'm looking for a way to visualize this data on a map.

Besides attaching a bar diagram to each link, I haven't found any way to achieve my goal using a GIS. The visualisation results should be printable so animations won't do.

#2 Example visualisation using diagrams:

All ideas welcome!


texts: how does projection happen in arcgis with regards to spatial reference

texts:
I have a map loaded in my Mapview called mMapview.and here is what i did
Point p1=new Point(66.9969,6.65428); //i took this coordinate from my gps device for demonstration
        Point p2=(Point) GeometryEngine
                .project(p1,
                        SpatialReference.create(4326),
                        mMapView.getSpatialReference());
        Graphic graphic = new Graphic(p2, new SimpleMarkerSymbol(Color.RED, 10, SimpleMarkerSymbol.STYLE.DIAMOND));
        mGraphicsLayer.addGraphic(graphic);//mGraphicsLayer is my graphic layer on my map view
when i run my program the marker was not there on my mMapview.but when i changed the statment  mMapView.getSpatialReference() with SpatialReference.create(32637) the marker appeared in exactly the right possesion i wanted it to be,so my Question is why this is happening?I checked my mMapview spatial reference id and latestid and and it is 102100and 3857 recpectively


texts: Is it possible to represent points as sparklines in ArcGIS for Desktop?

texts:
I have a point dataset with a number of attributes that I would like represent as sparklines (or simple line charts). This is similar to the bar chart symbology found within ArcGIS for Desktop. For example, if the points represent hydrocarbon wells and I have different pressure values in say 10 fields at different depths, the sparkline would show pressure fluctuations with depth across these values.

Does anyone know if this is possible?

The image shows a somewhat crude example.

enter image description here


texts: Visualizing elevation profile across a whole mosaic dataset?

texts:
I'm working with many DEMs derived from LiDAR data (IMG format) in a mosaic dataset. This allows me to create hillshade (or shaded relief) functions to display the data, which is great.

However it won't let me extract an elevation profile using the 3D analyst toolbar. Ideas?

I can get the profile tool to work on an individual DEM, but not across the mosaic dataset as a whole. (note: this dataset is too big to combine into a single raster).


texts: Counting (Multi)Polygons Within (Multi)Polygons on same Table

texts:
So I have a table with approximately 10k rows, with each row having its corresponding geometry entry. What I am attempting to do is count how many row geometries are contained/within the geometry of a single row as follows:

SELECT id, COUNT(*) FROM poly_table
WHERE WITHIN(poly_table.geom, poly_table.geom)
GROUP BY id ORDER BY COUNT(*) DESC;

However, all I get in return is a table with all rows having a COUNT(*) value of 1. So in understand they will certainly contain themselves but there are definitely polygons within polygons that are not accounted for.


texts: What is the metric of the radius PostGIS query?

texts:
I am trying to make a query (using raw query in sequelize orm for a PostgreSQL with PostGIS database) for a radius of a certain 'maxDistance' radius (preferably in meters):

const maxDistance = 100000000

const lat = 4.307662331310191
const lon = 52.09468958820517
const sequelize = hook.app.get('sequelizeClient')

return sequelize.query(`SELECT "id", "created_at", "wkb_geometry", ST_Distance(ST_SetSRID(ST_MakePoint(${parseFloat(lat)}, ${parseFloat(lon)}), 3857), "wkb_geometry") AS distance FROM "points" WHERE ST_Distance(ST_SetSRID(ST_MakePoint(${parseFloat(lat)}, ${parseFloat(lon)}), 3857), "wkb_geometry") < ${maxDistance}`)

When I create another point, i.e.: const HOME_LONLAT = [4.295980701608052, 52.09069818826458] in the database and run the query I get a result with a distance of 0.012344705180102 however the actual distance of these two geopoints is approx 930 meters. What am I doing wrong or how should I convert the distance to get the distance in meters?

EDIT:
When I alter the query into:

'SELECT "id", "created_at", "wkb_geometry", ST_Distance(ST_GeogFromText(\'SRID=3857;POINT(4.307662331310191 52.09468958820517)\'), "wkb_geometry") AS distance FROM "points" WHERE ST_Distance(ST_GeogFromText(\'SRID=3857;POINT(4.307662331310191 52.09468958820517)\'), "wkb_geometry") < 100000000'

I get this error:

error: Only lon/lat coordinate systems are supported in geography

EDIT2:
In PostGIS I have 3857 CRS:

SELECT ST_SRID("wkb_geometry") FROM points;
 st_srid
---------
    3857


texts: Reprojecting trouble: shapes saved in one EPSG load as another EPSG in QGIS and are 200m off?

texts:
I'm in a mixed working environment with some ESRI licenses and several people working with QGIS. For a project, we created data and received data in EPSG 31468 (Gauß-Krüger Zone 4). We now need to move to ETRS 1989 UTM 33N, seven digits.

I think EPSG 5650 is the one we need (ESRI users are told to use the apparently ESRI-specific WKID 102359).

I chose to use batch reprojecting the existing polygons, lines and points using ArcMap Batch Project, and produced output in EPSG 5650.

If I now create a QGIS 3.2.1-Bonn project with this EPSG, and load a created shapefile, it gets loaded as EPSG 5653.

From the linked description, I can see that only the following lines are different in EPSG 5650:

AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","5650"]]

I am no geographer, but I think the Easting/Northing definition is obsolete.
But why is it ignored? This seems to happen also with other projections.

The second part of my question is more severe:
comparing the GK4 data with, e.g., Google aerials shows they are correct. After reprojecting them, they are about 200m off. I suspect I need a transformation, but I don't know how I can do this. I can't choose one when reprojecting, neither in batch nor single mode.


texts:
I'm loading layer styles from .qml files in a Python plugin. I'm about to switch to loading them from a non-local PostgreSQL db. I've been looking at the (from QgsVectorLayer (v2.2))

listStylesInDatabase(QStringList &ids, QStringList &names, QStringList &descriptions, QString &msgError)
getStyleFromDatabase(QString styleId, QString &msgError)
loadNamedStyle(const QString &theURI, bool &theResultFlag, bool loadFromLocalDb)

For styles saved as .qml files loading them is straight forward. In terms of Python what are procedure for loading the styles from the db? Are they loaded simultaneously with the loading of the layer?
I have so far not found a Q&A on the forum regarding this.


texts:
I am trying to populate a field with the current time. Here is what I tried that did not work:

Start editing
Open Attribute Table
Add Field -> Name: time
    Type: Date
Right Click on Field Name (time)
    Open FieldCalculator
        Parser: Python
        Type: Date

Then I basically tried all of these expressions:

time.strftime('%d/%m/%Y')

That gives me the date, not the time

datetime.datetime.now( )

That gave me DATE and TIME. But as soon as I saved the changes, the time disappears.

(If the time would stay there I could do something with the FieldCalculator to populate another column like !time![:x])

datetime.timedelta(seconds=0)

That gave me <Null>.

Am I missing a step or something?


texts: Calculating real surface area using Zonal Statistics as Table Tool?

texts: I have a slope dataset.  I have a polygon shapefile with territory within the extent of the slope map.  how do I calculate the real surface area of this shapefile using the zonal statistics as table tool to get the sum of pixel values?  I am using ArcInfo 10.1


texts:
i want to make a 3d modell out of the mount everest and some other mountains. what i need therefore is a 14km² square around the peak from the mountain. so basically i use qgis to extract the area from the dowlnoaded SRTM geoTIFF (which i downloaded here http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp). qgis only needs the left upper and the right bottom corner of the areal, which i calculate with this site http://www.gpsvisualizer.com/calculators. So far, so good... now qgis cuts a rectangle!! but if i use microDEM to convert to a greyscale BMP i get a square shaped picture again!?!

the dem information of the bmp are the following:

Size: 13.92x13.95 km-------------> square (but not 14km²)

DEM size: 171x152----------------> no square

Screen Pixels: 517 x 518----------> almost square

Screen pixel size: 26.94 m

z range: 5040 to 8806 m

does anyone know what the heck is going on? i am really baffled and really dont know how to deal with this situation...


texts: OpenLayers Map Click Event not getting called in ReactJS

texts:
Problem: OpenLayers(v6.3) MAP is visible in REACTJS(16.13) but click event on MAP is not getting called.
Purpose: I want to place a marker and handle its click event to show some POPUP with data.
Code: Code included is though self explanatory; here is a quick explanation
[in ComponentWillUpdate()] MapComponent shown here receive {currentLocation} as property and moves map to that location. It then places a marker with custom Image. Its all working.
[in ComponentDidMount()] DIV ID reference is used to pass OLMAP reference.
I tried to register OLMAP click event in both ComponentDidMount() as well as ComponentWillUpdate() but it is not doing anything though it is getting executed as last console message gets printed.
Only problem in this code is OL-MAP CLICK event Handling. What is wrong?
    import React, { Component } from 'react';
    import ReactDOM from 'react-dom';
    import Map from 'ol/Map';
    import View from 'ol/View';
    import TileLayer from 'ol/layer/Tile';
    import OSM from 'ol/source/OSM';
    import VectorLayer from 'ol/layer/Vector';
    import VectorSource from 'ol/source/Vector';
    import Feature from 'ol/Feature';
    import {Style,Icon} from 'ol/style';
    import Point from 'ol/geom/Point';

    import { fromLonLat } from 'ol/proj.js';
    import LOCPNG from '../../../assets/images/location.png';
    class MapComponent extends Component {

      constructor(props) {
        // console.log('[MC.js] constructor');

    super(props);
    this.olMap = null;
    this.mapRef = null;
    this.setMapRef = element => {
      this.mapRef = element;
    }
  }
  componentDidUpdate(prevProps, prevState) {
    console.log('[MC.js] CDUpdate');
  }
  componentWillUpdate(nextProps) {
    // console.log('[MC.js] CWU');
    const mapView = this.olMap.getView();
    mapView.animate({
      center: nextProps.currentLocation,
      duration: 2000,
      zoom: 10
    });

    this.makeMarker(nextProps.currentLocation, 0);

    this.olMap.on('click',(evt) => {
      console.log(evt);
      console.log(evt.pixel);
      const feature = this.olMap.forEachFeatureAtPixel(evt.pixel, (feature) => feature);
      if(feature){
        console.log('PopUp');
      }
    });
  }
  makeMarker = (location, index) => {
    let marker = new Feature({
      geometry: new Point(location)
    });
    marker.setStyle(new Style({
      image: new Icon(({
        // crossOrigin: 'anonymous',
        src: LOCPNG,
        // enter code here`enter code here`
      }))
    }));
    let vectorSource = new VectorSource({ features: [marker] })
    var markerVectorLayer = new VectorLayer({
      source: vectorSource,
    });
    this.olMap.addLayer(markerVectorLayer);
  }

  componentDidMount() {
    // console.log('[MC.js] CDM');
    const mapDOMNode = ReactDOM.findDOMNode(this.mapRef);
    // console.log('[MC.js] CDM',mapDOMNode);
    this.olMap = new Map({
      target: mapDOMNode,
      layers: [
        new TileLayer({
          source: new OSM()
        })
      ],
      view: new View({
        center: this.props.currentLocation,
        zoom: 10
      })
    });

    this.olMap.on('click',evt=>{
      console.log(evt);
      console.log(evt.pixel);
      const feature = this.olMap.forEachFeatureAtPixel(evt.pixel, (feature) => feature);
      if(feature){
        console.log('PopUp');
      }
    });

    console.log('[MC.js] CDM');
  }

  render() {
    // console.log('[MC.js] render()');
    const styles = { height: '100%', width: '100%' }
    return (
      <React.Fragment>

      </React.Fragment>
    )
  }
}

export default MapComponent;


texts: Which Tools do you use for Classification of Remote Sensing Data?

texts:
Which tool do you prefer to use for classification of remote sensing data, e.g. classifying land use, and why?

Which other tools have you tried, and why did you decide against them?


texts: Is it possible to use the Goode Homolosine projection in QGIS? Both QGIS and GRASS seem to have trouble with the interruptions, squashing the higher latitudes together. I tried +proj=goode +ellps=sphere +datum=wgs84 as a custom CRS.


texts: How to copy the location of these markers on a website Google map?

texts:
I want to know the location coordinates of all those points in order to put them in GIS software. (QGIS)
Is there an easy way to do it?
This is the website: https://www.ah.nl/suppliers/map


texts: How do you filter a layer on a joined table in QGIS?

texts:
I am attempting to do a join in QGIS between a SQL Server Spatial table and a non-spatial table (on a single integer field).  I was able to add both successfully and make the join, but the query builder dialog only shows the fields from the spatial table.

I also tried setting a feature filter on the non-spatial table before joining to the spatial table, but that filter was not honored (no error, but no affect on the layer).

QGIS version: 2.12.1-Lyon    QGIS code revision: 16760fd

Microsoft SQL Server 2012 - 11.0.2100.60 (X64)
    Feb 10 2012 19:39:15
    Copyright (c) Microsoft Corporation
    Enterprise Edition (64-bit) on Windows NT 6.1 <X64> (Build 7601: Service Pack 1) (Hypervisor)


texts:
I am doing Project on Drainage basin, to find out the Elongation ratio I need to find out Basin length, Can any one of you tell me how to find out BASIN LENGTH using ArcGIS or any other way. My data base is like below Image

(i just uploaded this image from internet)


texts:
I want to share data on a project with many users, and wanted something similar to ArcReader. I like this software as it can be accessed from the server for speed, and therefore updated to other offices periodically.  If there was an option to add layers for users to mark up that would be awesome.

Is there a similar solution or approach for QGIS?

I realize that the users could just use QGIS but was hoping for something so straight forward it could be used virtually straight away.


texts: In GEE how to select bands for analysis base on availability

texts:
I'm trying to use a simple threshold to detect water in Sentinel SAR data.  For my application, the VH polarization is more reliable, but I can use the VV polarization if VH is not available. What's the best way to write code to use the VV band if the VH is not available?
I've created a simplified version of my code here to illustrate what I'm trying to accomplish.  It involves the problem of mixing server side data with client side processing and using an 'if' statement. So,how should something like this be accomplished?
var genWaterMask=function(image,threshold){
    return image.lt(threshold);
  }

var aoi=ee.Geometry.Point(-105.68, 28.12).buffer(9000);
var boundsSet = ee.ImageCollection('COPERNICUS/S1_GRD')
                .filter(ee.Filter.eq('instrumentMode', 'IW')).filterBounds(aoi);

var fall = ee.Filter.date('2015-09-21', '2015-12-20');

var bandName='VH';
var threshold=ee.Number(-10.0);

var fallSet=boundsSet.filter(ee.Filter.listContains(
  'transmitterReceiverPolarisation', bandName)).filter(fall);

if ( fallSet.size() < 1 ){
  var bandName='VV';
  var threshold=ee.Number(-15.0);
  var fallSet=boundsSet.filter(ee.Filter.listContains(
    'transmitterReceiverPolarisation', bandName)).filter(fall);
}

var waterMask=genWaterMask(fallSet.select(bandName).mean(),threshold);


texts: How to use more CPU while processing in QGIS

texts:
I'm trying to do some clip process with a quile large layers but the process is too slow. When I see the task administrator, just around 28% of the CPU is being used.
I've tried with 32-bit and 64-bit versions of QGIS 2.8 and 2.12. Same results: A clip of this layer takes hours, while ArcGIS 10.2 just needed 2 minutes.
Does someone know what could be happening ??


texts: Is there a way to put ArcMap's attribute table on a different thread or processor?


texts: Converting SpatialGridDataFrame into polygon shapefile in R?

texts:
I have a SpatialGridDataFrame in R and want to plot just one of its attributes as as polygons in a GIS. There seems to be no simple way to do this.

The situation is complicated by the fact that the variable I want to export is not numeric so I can't just convert it to a raster and export it. It seems like I have to have the grid file converted to polygons in R and then export the polygons into a GIS format.

I've tried Grid2Polygons package but I don't seem to be able to carry across the attributes.


texts: Efficient way to check sample points on many Landsat scenes?

texts:
I want to do an accuracy assessment for an annual burned area map (large area - whole country). I have random sampling points which I need to check in simple burned - not burned categories but throughout the fire season of several months. I have sample points in around 60 Landsat scenes and depending on cloud cover there will hopefully be imagery available from several acquisition dates. Assuming that it will be possible to use 3-4 dates per scene, I need in total at least 180-240 scenes.

Is there an efficient way/online tool/software to look at Landsat-8 scenes in different band combinations without having to download them? Something where I can use my point layer and click through an archive of Landsat data?

I fear I am dreaming of something non-existent...any hints?


texts: Changing units in QGIS: degrees to meters

texts:
I have a global map with CRS WGS84-EPSG:4326 and I want to calculate the area of my polygons. From what I read, to change the unit from degrees to meters I need to change my CRS. However, there is no global system for this.
Is there any solution to change the unit from degrees to meters for a global map?


texts: Error while adding WMS layer in QGis 2.18.13

texts:
I already set up a geoserver instance and add some raster stores to it. I want to use them as WMS service through different aplication like QGis and Autocad but when i try to add the WMS layer in QGis 2.18.13 i get this error:

Invalid Layer: WMS provider Cannot calculate extent Raster layer Provider is not valid (provider: wms, URI: IgnoreAxisOrientation=1&IgnoreGetFeatureInfoUrl=1&IgnoreGetMapUrl=1&contextualWMSLegend=0&crs=EPSG:32634&dpiMode=7&featureCount=10&format=image/png&layers=raster_utm:layer1_utm2&styles=&url=http://localhost:8080/geoserver/ows?version%3D1.3.0%26

This does not hapen in Qgis 3.0 or in Autocad map 2017. The EPSG used for the layer is 32634.


texts: Using ModelBuilder to select layers by attributes, dissolve and then merge them?

texts:
I would like to reclassify polygons, I want to combine smaller polygons into bigger polygons based on their ID (each polygon has a unique ID). I want to regroup ~800 polygons into ~60 polygons.
So in the first round I, e.g, want to select polygons with the ID 3,5,7 and dissolve them into one. In the second round I want to select polygons with the ID 34,666, 45,6,77 and dissolve them...until all the polygons are regrouped and dissolved. Then I would like to merge all the dissolved polygons into a new shapefile with only ~60 polygons.

So I thought I can automate this with the ModelBuilder, first I tried an approach with

Make Feature Layer then Select by Attributes then Dissolve.
But that way I still have to change the SQL statement in the Select Layer By Attribute each time.
Also I wouldn't know how to save the dissolved output layer after each run so that it doesn't overwrite it. Is there a function that keeps merging the previous model output with the current output?

While searching for solutions online I came across the Feature Class to Feature Class option as it apparently makes the Make Feature Layer and Select Layer By Attribute steps unnecessary. How could I then use a list variable of SQL statements and a list variable of string names for the different outputs as input?


texts: Leaflet - Adding another layergroup on top of a heatmap layer

texts:
Layering with z-index works when using layergroups - but is there any working way to have another layer with circles or other elements ontop of a layer with a heat-map?

The heat-map seems to be always on top - except for single markers (circle-markers also don't work).

Using Mapbox/leaflet and https://github.com/Leaflet/Leaflet.heat


texts: Wrong reprojecting of a shapefile with to_crs in GeoPandas

texts:
When reprojecting an ESRI Shapefile (EPSG:21781) into an EPSG:4326 with QGIS and then importing the resulting file in GeoPandas and plotting it, the coordinates are correct.

When importing the same shapefile directly in GeoPandas and then reprojecting it using to_crs({'init', 'EPSG:4326'}), there is a slight offset.

Here is the result of my_shapefile.crs:

{
  'ellps': 'bessel',
  'k_0': 1,
  'lat_0': 46.95240555555556,
  'lon_0': 7.439583333333333,
  'no_defs': True,
  'proj': 'somerc',
  'units': 'm',
  'x_0': 600000,
  'y_0': 200000
}

I don't understand if I missed a step with the CRS conversion in GeoPandas or if the to_crs function needs another parameter.


texts: Python Geolocation Software to Filter Cities By Population Within a Radius?

texts:
What I'm trying to achieve

I want to enter a US city/state or zip code value.  Then enter a radius (ex: 30 - would stand for 30 miles).  Then I enter a minimum population (ex: 100000 - would tell the program logic to only look for cities with populations greater than or equal to 100000).  The program returns the names of cities within the radius that has populations greater than or equal to the minimum population.

I don't care about super-precise numbers.  For example, if a major city barely touches the outside radius, it still counts so I want to include that city's entire population (even though it's only perhaps a small sliver of the actual population) in the output.

My Question

Is there an open-source Python software or solution out there that would allow me to achieve this?  If yes, what do you recommend?

What I'm considering

Google geolocation but don't know if it's possible to do what I'm trying to do
Some Python radius software package I don't know about
Utilizing a giant database lookup table or an API but don't know which route is best from a strategic standpoint?


texts:
I am a new Arc Objects user, with little ESRI experience.

We are using VBA to create a spatial report on a Map Control using a mdb table. The mdb table contains a Lat and Long field, and this data is collected in the Access application.

For the reporting, we need to display a map with the locations from the mdb table.

My questions are:

How do I access data in a mdb table, as spatial data? Does it have to be registered as spatial? I am not sure if this can be done at runtime or has to be done using Arc Catalog.
At runtime, if the data from the mdb table is in a spatially recognized table, do I just create a new layer, set the symbology, and then set a "where" clause on the new layer?


texts: Set PATH variables for ogr2ogr without admin privileges

texts:
In the office I use ogr2ogr to convert several file formats such as KML or ESRI shapefiles. Since the commands are the same almost every time I'd like to run ogr2ogr with a batch. I'm using the GDAL package for Windows from gisinternals. Normally running the SDKShell.bat would set all the PATH variables that are needed to use ogr2ogr through cmd but our IT-department blocks users from changing the PATH. Therefore I can't call ogr2ogr with a batch because cmd doesn't know where to look for it.
My current approach is to start the SDKShell.bat and copy my commands into it each time. But that is not really what i wanted to achieve. Is there a possibility to solve this issue?


texts: Is there a way to define the relative position of items in print composer?

texts:
Is there any chance to make the items in print composer relative to the margins. Let's say for example the legend. I normally align it to the corner on the right bottom where the legend is 10 px away from the boarders of the map.

If I add a new item to the legend, the legend gets bigger and I have to align it again. To give the legend a relative position independently from the legends size would be very favourable (tell qgis to always align the legend in the corner on the right bottom with a distance of 10 px always). Than I could define one layout and make several maps without readjusting those items all over again...


texts: Do the Data Interoperability Tools need a Special Licence?

texts:
I am tring to add a function quick export from Data Interoperability Tools to a Model Builder in ArcGIS10.0 but I have an error called Error 000816: the tool is not valid.

That is because I need a license or what is the problem and how can I solve it?


texts: GetFeatureInfo for multiples layers openlayers 6

texts:
I'm trying to get information from several layers in Openlayers 6, but it only returns information from one. I define the layers and a source with the url and the layers as params
var spDo = new ol.layer.Image({
    title:'DO's',
    baseLayer:false,
    source: new ol.source.ImageWMS({
       url: 'http://localhost:8080/geoserver/grap/wms',
       params: {'LAYERS': 'land:do_s'},
       serverType: 'geoserver',
    }),
    visible:true,
    active:false,
    //opacity: 0.4,
    displayInLayerSwitcher:false
});
var spParc = new ol.layer.Image({
    title:'Parcels',
    baseLayer:false,
    source: new ol.source.ImageWMS({
        url: 'http://localhost:8080/geoserver/grap/wms',
        params: {'LAYERS': 'land:geo_land'},
        serverType: 'geoserver',
    }),
    visible:true,
    active:false,
    //opacity: 0.4,
    displayInLayerSwitcher:false
});
var spParCom = new ol.layer.Image({
    title:'Control Parcels',
    baseLayer:false,
    source: new ol.source.ImageWMS({
        url: 'http://localhost:8080/geoserver/grap/wms',
        params: {'LAYERS': 'land:vw_parcels'},
        serverType: 'geoserver',
    }),
    visible:true,
    active:false,
    //opacity: 0.4,
    displayInLayerSwitcher:false
});
var demRequest= new ol.source.ImageWMS({
    url: 'http://localhost:8080/geoserver/grap/wms',
    params: {'LAYERS': 'land:do_s, land:geo_land, land:vw_parcels'},
});
map.on('singleclick', function (evt) {
    document.getElementById('info').innerHTML = '';
    var viewResolution = /** @type {number} */ (view.getResolution());
    var url = demRequest.getFeatureInfoUrl(
        evt.coordinate,
        viewResolution,
        EPSG:3857',
       {
        'INFO_FORMAT': 'text/html',
        'QUERY_LAYERS': 'land:do_s, land:geo_land, land:vw_parcels
       }
    );
    if (url) {
        fetch(url)
         .then(function (response) { return response.text(); })
         .then(function (html) {
        document.getElementById('info').innerHTML = html;
    });
    }
});

I don't understand the problem, I define the list of the layers in the url params and in the 'QUERY_LAYERS' according to the documentation.
The server response contains the three layers:
 ...REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=land%3Ado_s%2Cland%3Ageo_land%2Cland%3Avw_parcels&LAYERS=land%3Ado_s%2Cland%3Ageo_land%2Cland%3Avw_parcels&INFO_FORMAT=text%2Fhtml...


texts: Loading online OSM raster tiles using URL in custom QGIS application (C++)

texts:
I am trying to write a C++ version of a custom QGIS (v 3.10) application I initially wrote in Python that overlays data onto OpenStreetMap.  I am able to do this successfully by loading online tiles via URL as a raster layer in Python using WMS provider; however the same URL does not work when using the C++ API, despite being able to load other raster files (e.g. GeoTIFF).
Code in Python:
url = 'type=xyz&url=https://a.tile.openstreetmap.org'
url += '/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0&crs=EPSG3857'
prj = QgsProject()
qmc = QgsMapCanvas()
layers = []
ras = QgsRasterLayer(url,'OpenStreetMap','wms')

if ras.isValid():
    print("Basemap loaded successfully!")
    prj.instance().addMapLayer(ras)
    qmc.setExtent(ras.extent())
    layers.append(ras)
    qmc.setLayers(layers)

 else:
     print("Unable to load basemap.")

C++ version:
QString url = "type=xyz&url=https://a.tile.openstreetmap.org";
url.append("/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0&crs=EPSG3857");
QgsProject() *prj = new QgsProject();
QgsMapCanvas *qmc = new QgsMapCanvas();
QList <QgsMapLayer *> layers;
QgsRasterLayer *ras = new QgsRasterLayer(url,'OpenStreetMap','wms');

if ( ras.isValid() )
{
    qDebug() << "Basemap loaded successfully!";
    prj->instance()->addMapLayer(ras);
    qmc->setExtent(ras->extent());
    layers.append(ras);
    qmc->setLayers(layers);
 } else
{
     qDebug() << "Unable to load basemap.";
 }

I get the successful message and loaded map in Python and the invalid message and no map in C++. I have tried using other URL's but nothing has worked thus far.  If there are plugins available that might help, but I would still need to write workable code without the use of the QGIS gui itself (as in most examples using plugins).
EDIT Using the advice of @ThomasG77 below, I made the following amendment:
QStringList plist = QgsProviderRegistry::instance()->providerList();
for ( int i=0; i<plist.length(); i++ )
{
    cout << plist.at(i).toStdString() << endl;
 }

which yields the following output:
gdal; memory; mesh_memory; ogr

Do I have to manually add the OSM/WSM provider in my .pro file or environmental variable?  Using Fedora Linux (31)


texts: Map Matching (snap GPS points road) in ArcGIS (REST API)?

texts:
Does ArcGIS have Map Matching function? (In Google Maps Services it's called "Snap to roads" - https://developers.google.com/maps/documentation/roads/snap)

I checked next link, but didn't find anything: https://developers.arcgis.com/rest/services-reference/get-started-with-the-services-directory.htm

What is Map Matching: https://en.wikipedia.org/wiki/Map_matching

Example:


texts: Using legal descriptions of lands and their subsections to draw polygon?

texts:
I have Legal descriptions of land (Section, Township, Range) and their subsections (N1/2, S1/2, E1/2, W1/2, NW1/4, SW1/4, N1/2NW1/4, etc.)

I want to have a polygon encompass the subsections of the legal plots.

I intend to use county maps or a state map with already-built-in legal descriptions, so I believe there shouldn't be an issue with scale - the legal description should take me to the exact plot I want to draw over.

I'm working with a high volume of data so ideally I'll automate the drawing of the polygons.


texts: Filling holes by acres threshold using PostGIS?

texts:
I have a layer that has holes in it and I would like to fill only the holes that are smaller than 1/10th of an acre

I know how to fill all the holes with this code

SELECT id, ST_Collect(ST_MakePolygon(geom)) As geom
FROM (
    SELECT gid, ST_ExteriorRing((ST_Dump(geom)).geom) As geom
    FROM layer
    ) s
GROUP BY id

but how I was not able to figure out the acreage of the holes.

should I be utilizing https://postgis.net/docs/ST_DumpRings.html ?

not sure how to think this process through.

I was thinking use the above query to fill in all the holes than run an st_difference from the swiss cheese hole layer to the filled holed layer and the output would be those holes as polygons where I can calculate the acreage..but then what merge them back to the swiss cheese hole layer?

is there a simpler approach?

I ended up running this query

select geom,st_area(geom)/43560 acres from(
    select (st_dumprings((st_dump(shape)).geom)).geom geom
    from agriculture.agricultural_grtr_250) t
where st_area(geom)/43560 > .1

but it did not yield the results I wanted. I wanted to keep the holes greater than .1 but it still filled them.

I guess I want the SQL code that will emulate the ArcGIS tool eliminate polygon part


texts:
I'm running 2.14.11 and I'm getting the following error when trying to install TimeManager

Plugin installation failed:Failed to unzip the plugin package. Probably it's broken or missing from the repository. You may also want to make sure that you have write permission to the plugin directory:

Is this plugin no longer available?


texts: Model builder - Create a new File GDB then convert a CSV into a file GDB table in new File GDB

texts: I am trying to create a model to create a new file GDB.  The model would then iterate through a folder of CSV tables, convert them to file GDB tables and add them to the newly created File GDB.  I need to link the newly created file GDB to the copy rows function so it automatically puts any tables into the new file GDB.


texts: QGIS Processing modeler not showing toolbar

texts: I've just installed QGIS recently (3.10.10 with GRASS 7.8.3) on a new computer. Previously, I used the same version and there's a toolbar in the Processing Modeler (to save and run model). But now, I cannot see it anymore. Why is it?


texts: Adding Spatialite table with clip results in QGIS?

texts:
I have 2 layers and I want to make a clip of layer 1 with layer 2 in SpatiaLite (I'm using the DB Manager within QGIS).

That's the code I'm using:

create table clipped_lines as
select * from lines, estensione
where st_intersects(lines.GEOMETRY, extension.geom)

The table is created but it is not a geographical object so I cannot add it to the QGIS map canvas.

The problem (I think) is that the geom column is seen as a normal column and not as a GEOMETRY column. I've tried to update the column with:

ALTER TABLE clipped_lines
ALTER COLUMN geom TYPE GEOMETRY(LineString,4326);

but nothing happened.

I thought that the st_intersects function updates automatically the geom column.


texts: Drawing points in real time (or at short intervals) using PyQGIS

texts:
I’m working on a robot that moves in space and transfers its GPS coordinates to the server. I need to create a script in QGIS that will receive coordinates and draw points based on these coordinates in real time. I have a script that draws points from a list and I know how to make a client that will receive coordinates from the server. But I don’t know how to make a script that will work constantly and in real time draw the received coordinates (or work cyclically and draw points at a small interval, this solution will also work).
How can I solve this problem?


texts: How to edit the capabilities with python when publishing map services

texts:
Publishing map service the normal way is done here using arcgis desktop.

But I decided to do it with python so I can automate the process or something like that but first off I need to edit some of the info according to this, you would create the SD file before publishing.
Due to the 0 rep I can't post more then two links but normal publishing

and using arcpy

.

As you can see I highlighted the description according to the image does it means it is searching for the Item Description.

Is the the right way of doing it

# find the Item Information Description element
Capability= doc.getElementsByTagName('Capabilities')
for item in Capability:
    if item.parentNode.tagName == 'FeatureAccess':
       #How do i check this check box?
    else:


texts:
I need to create a raster with values showing directions relative to a single point. Each cell will have the azimuth to the origin point.

Basic trigonometry could calculate this using the differences in x and y coordinates of the 'point.shp' for each cell but I am not aware of a suitable QGIS process.

Is it possible to access the geographic coordinates for each cell within raster calculator, to allow me to calculate the azimuths for each cell into a new raster?


texts:
I want the title (a label) of my map in print composer to automatically contain the name I've given the print composer and also the date that it was created. For example I want it to say 'Really cool map of (print composer name) created on (date the print composer was created).'

I can then use it as a template for lots of other maps which will automatically update with the site name and date created.

I think I need to use the 'Insert an Expression' but I can't figure it out.


texts:
Im running Qgis server on debian. It’s a virtual machine.

i already tried searching, but i found nothing, what helped me out.
What i did:
I installed QGIS Server and all needed packages. I can talk to the server via http, when i call this link:
http://IP/cgi-bin/qgis_mapserv.fcgi

Answer:

   <ServiceExceptionReport version="1.3.0">
<ServiceException code="Service configuration error">Service unknown or unsupported</ServiceException>
</ServiceExceptionReport>

So far, so nice.

When i call:

http://IP/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&MAP=/var/www/html/S5/qgis/test/qgis-demo/world3.qgs

i get:

    <ServiceExceptionReport version="1.3.0">
<ServiceException code="WMS configuration error">
There was an error reading the project file or the SLD configuration
</ServiceException>
</ServiceExceptionReport>

I set the rights to the Folder /var/www/html/S5/qgis and all Folders and files inside to 775 via chmod but i cant get it to work. I have no idea and i find no usefull threads in the hole internet to help me out. Problem 2 is, that qgis-server is not writing any logs.
For logs, i set up a config in

/etc/apache2/mods-enabled/fcgid.conf

like this:

#Tell QGIS Server instances to use a specific display number for xvfb
# necessary for printing, see below
FcgidInitialEnv DISPLAY ":99"
# Activate QGIS log (different from apache logs)
FcgidInitialEnv QGIS_SERVER_LOG_FILE /var/log/qgis/qgisserver.log
FcgidInitialEnv QGIS_SERVER_LOG_LEVEL "0"
FcgidInitialEnv QGIS_DEBUG 1
# Add a default QGIS project
#SetEnv QGIS_PROJECT_FILE /var/www/html/S5/qgis/qgis-demo.qgz
# QGIS_AUTH_DB_DIR_PATH must lead to a directory writeable by www-data
FcgidInitialEnv QGIS_AUTH_DB_DIR_PATH "/var/www/qgis-server/qgisserverdb/"
FcgidInitialEnv QGIS_AUTH_PASSWORD_FILE "/var/www/qgis-server/qgisserverdb/qgis-auth.db"

<IfModule mod_fcgid.c>
# Longer timeout for WPS... default = 40
FcgidIOTimeout 120
FcgidMaxRequestLen 26214400
FcgidConnectTimeout 60
</IfModule

The owner of the logfile is www-data.
I set a # to comment the standard projct file out to serve the WMS via the attribute MAP.
I tried other project as well, but its not working.
Can anybody help me out please?


texts: QGIS DB Manager import Primary Key issue

texts:
I have a number of layers that I want to import into a GeoPackage using DB Manager > Import Layer.  This data will end up in a process that controls data entry and it needs to comply with corporate naming conventions - needs to be "gid".  On the import screen I tick the Primary Key and entered "gid".  The import tool successfully imports the data, but ignores anything I enter as the PK, it uses "fid" every time.

I have tried different options, "gid", "kid", even left this unticked - same result = "fid".
None of my tables have "gid" as an existing attribute field.

Is this a bug and is there a workaround? BTW, I use QGIS 3.10


texts: How to find the nearest point projected on the road network?

texts: I would like to know if there is a function at pgRouting that given a point (lon/lat) can give you the projection of that point to the road network .


texts: Why is the mxd property "Pathnames" disabled?

texts:
I have an .mxd that is using absolute pathnames. I want to change this property to use relative pathnames but the option is disabled, grayed out, when I go to the .mxd properties. It is actually grayed out every every .mxd I try. How do I enable this option?

Edit

At @kenbuja suggestion I checked the option "Make relative paths the default for new map documents" in Arcmap options. Then I was able to get new .mxds to use relative paths but the option was still disabled in the mxd properties. Also I need to try to get this working about 20 existing .mxds.


texts: Using numpy.dstack() for merging two 3D arrays but losing spatial reference.

texts:
The code is as following:

>>> temp = numpy.dstack((spi_negative,spi_positive))
>>>
>>> temp.shape
(16, 126, 294)
>>> out_ras = arcpy.NumPyArrayToRaster(temp,pnt,cellsize1,cellsize2,value_to_nodata=0)
... outfile = "D:/TRMM/Gurminder_SPI_TOTAL.img"
... out_ras.save(outfile)
... arcpy.DefineProjection_management(outfile, spatialref)

The shape is maintained and raster image is saved successfully. But the result shows both the stacked images with different position.


texts: Can I have a layer in my map that gets its imagery from say google? (ArcGIS Mobile 10)

texts: I've been messing with ArcGIS Mobile and I get that it loads and caches layers from an ArcGIS Server for disconnected use. I'm wondering if it supports layers where data is loaded from a layer in realtime? I can't find any examples where I create a layer, give it a URL and have it load imagery data as I pan the map.


texts: Define Projections for layers with multiple locations in MXD's

texts:
# Script to use Layers in MXD source file paths to locate and redefine projections

import arcpy, glob, os, sys, arcgisscripting

mxd = arcpy.mapping.MapDocument("N:\MXDs\_ThomaKL\MXD GPS Valves\GPS_Shapefile_Paths.mxd")

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        print lyr.dataSource

from arcpy import env, mapping
path = lyr.dataSource
env.workspace = path
env.overwriteOutput = True

print 'Reading files from ' + path

sr = arcpy.SpatialReference(103133)

x=0
z=x+1
NoProjCount=0
FileList= arcpy.mapping.ListLayers(lyr)
for lyr in FileList:
    desc = arcpy.Describe(lyr)
    SR = desc.spatialReference
    if SR.name == "Unknown":
        print "Projection of " + str(lyr) + " is " + SR.name + " so defining projection."
        f = open('NoProjection.txt', 'a')
        f.write(str(lyr)+"\n")
        f.close()
        arcpy.DefineProjection_management(lyr, sr)
        NoProjCount=NoProjCount+1
    else:
        print lyr + " is projected " + str(SR.name)

so I have a script that works, but it only defines the projection for the very last Layer in the MXD and not the entire list of layers (there are around 50). I can tell I am missing something but I am not sure what. below is what the NoProjection.txt file has stored in it.

  Valves.shp GPS Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points GPS
  Points\DataCollectedLater-Valves\Valve_GPS Points


texts: Splitting polygon into smaller rectangular polygons that have user input dimensions (length x width) using ArcPy?

texts:
Given any irregularly shaped large polygon, I'd like to be able to subdivide it into as many (full) smaller rectangular polygons as possible (like a grid, but rectangular and with user-specified dimensions).

I know of the Tesselate function in ArcGIS Pro, and of the Fishnet function, but I'm not totally sure those would have the outcome I'm after.

I can manage writing something in ArcPy to facilitate this, I'm just not familiar with the built-in functions in Arc well enough to know what functions to look into leveraging.

The workflow I'm envisioning is something along the lines of:

Create a grid that has the user-specified cell dimensions (a x b)
Clip the grid to the size of the polygon, and create shapefiles for each cell.

All of the space from the original polygon needs to be accounted for in the final output.

A good example use case for this would be, you have an area of irregularly shaped  parcel of land and want to subdivide it into equal rectangular lots and create polygons for each of them.

Have access to ArcGIS Pro, ArcGIS 10.3, and can build arcpy toolboxes if needed.

Image:
 - "Grid" in red
 - Original irregularly shaped polygon in black.
 - Output polygons in blue


texts:
I'm currently trying to do a connectivity analysis of roads of a specific type. (how far can you get before encountering a dirt road?) Ideally, I want the map to look something like this:

It is basically a connectivity analysis by field. So each "web" should be merged as a single feature because the roads layer I have now is a bunch of independent road segments(features). I would do it manually merging all the features, but there are over 5000 features in this roads layer.
Is there any tool that can help?


texts:
I have line layer as shown in image 1. I created middle point of line by using the below query in virtual layer.
SELECT ST_StartPoint(ST_Line_Substring(geometry,0.5,1))
FROM "TR_ROAD copy"

The problem is it only creates the middle point without creating data (ids) for middlepoints(nodes).

Is there are way to create unique ids for each middle point as sample is shown in last image


texts: Performing operations on Python Script tool startup

texts:
I want to perform a selection on tool startup in order to obtain a choice list for parameter derived from selected features. Code kind of works, as I get my choice list with correct values, however it works in a strange way:

After pasting code into validation window and clicking OK, selection starts and is performed several times.
On the first tool startup, selection is done again several times.
Choice list gets remembered as Value List and doesn't change with further tool startups (in fact on next startups, nothing new happens - no selection or any other actions from initializeParameters).

What am I doing wrong?
Here is a code:
import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

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

def initializeParameters(self):
  """Refine the properties of a tool's parameters.  This method is
  called when the tool is opened."""
  arcpy.SelectLayerByLocation_management(in_layer="OBR",   overlap_type="CONTAINS", select_features="POW")
  choice_list = []
  with arcpy.da.SearchCursor("OBR", ["OBR", "NAZWA"]) as cursor:
      for row in cursor:
          choice_list.append(row[0])
  arcpy.SelectLayerByAttribute_management("OBR", "CLEAR_SELECTION")
  self.params[0].filter.list = choice_list

  return

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

  return

def updateMessages(self):
  """Modify the messages created by internal validation for each tool
  parameter.  This method is called after internal validation."""
  return


texts:
For more than two days I'm trying to make the TSP algorithm work in PostgreSQL/PostGIS. My problem seems to occur when I'm trying to take into account the road network for TSP algorithm. I know  that the algorithm it will calculate the distance in Euclidean space as a default setting.

As an example:

SELECT vertex_id
FROM tsp('SELECT id as source_id, lon AS x, lat AS y
          FROM point_layer','23,221,63,432,335,',63);

where the point_layer is the table which holds the points information (as PlaceMarks), numbers like: 23,221,63,432,335 stands for places (cities) and 63 as a starting point.

If you have experience with this approach please advise.


texts: How to convert a spatialpolygon to a SpatialPolygonsDataFrame and add a column to the attribute table

texts:
coast<-readShapeSpatial("coastline.shp")
landc<-readShapeSpatial("landcover.shp")
ro<-readShapeSpatial("roads.shp")
bc<-gBuffer(ro,width=100)
landc$ratings=1
landc$ratings[landc$LANDUSE_ID==4]=0

Above, I take any category that has 4 and in the new column put it as 0.

At this point, I want the column named ratings as well for the bc, where it will take 0 if it is inside the buffer and 1 if it is outside. The problem is that the bc is SpatialPolygons and it doesn't contain the attribute table.

Obviously to add a column to a SpatialPolygon object you have to convert it to a SpatialPolygonsDataFrame, but I don't know how.

I tried this:

buf_df<-as.data.frame(bc)
s_po<-SpatialPolygonsDataFrame(bc,buf_df)
s_po$ratings=0

but this error pops up:

row.names of data and Polygons IDs do not match


texts: Publishing Hosted Feature Layers to ArcGIS Online?

texts:
If you publish hosted feature layers, is there an option to make them private afterwards?

Reason being I have a layer I want to use using ArcGIS Online but it says the file is too large and the other option is to publish it. If there are any other workarounds I'd be happy to know.


texts: I'm new to arcgis javscript api, and trying to get layers from a group layer. I'm not sure how to get this done. any help?


texts: Excluding layers from snapping using IPointSnapper

texts:
I'm working with ArcObjects in C#, and am trying to exclude all map layers but one (the only one I want to snap to) from the snapping environment. I'm using the ISnappingEnvironment and IPointSnapper interfaces. I've put the following code in the OnClick method of the tool I'm using:

 m_Snapper = m_SnappingEnvironment.PointSnapper;
 ISet anExclLyrSet = new SetClass();
 m_Snapper.ExcludedLayers(ref anExclLyrSet);

 // ...loop through all map layers, adding all but the one of interest to anExclLyrSet

 m_Snapper.ClearCache();

Basically, I'm getting poor snapping performance due to the presence of a basemap layer that contains multiple feature layers. If I turn that layer off (make invisible), the performance is good.

I've confirmed that all the layers (for composite type layers, I add each sub-layer to the exclusion as well) are being added to the exclusion set, but I can't seem to get the performance to improve without turning off that basemap layer. I'm not sure what I'm doing wrong - does the cache need to be cleared after each individual layer is excluded? I'd appreciate insights anyone might have regarding this. Thanks!

Mark C.


texts: I am trying to import a historic USGS map illustrating an island in the Mississippi River from 1953. I would like to project the scanned map onto a basemap without georefrencing with control points but with the lat/long coordinates given on the map. However, I cannot seem to figure out which coordinate system to use. I have tried NAD27 and many other, but when i display it on a basemap, it ends up being somewhere in the middle of the atlantic or Caribbean. Attached is the information given on the map.


texts: Digitizing polygons against georeferenced map using ArcMap

texts:
I am drawing polygons of counties over a 17th Century map that I have georeferenced. I have an issue with this version of ArcMap which I haven't experienced with past versions.
The shape keeps closing after 3 vertexes have been applied.  The only way I have found to get round this is to draw a square of vertexes over the whole county and then keep adding vertexes to draw in the polygon which is very time consuming.
How can I avoid this happening?


texts: Modify a script in the R Script repository in QGIS

texts:
I'm trying to use to R scripts available in QGIS via the processing Panel, and in the "Get R scripts from on-line script collection". I spotted some errors in the script "Kriging" that prevents it from working. I corrected the errors to make it work, but I think it would be nice to correct the source script with the correction I made.

Does anyone knows how can I update a script from this list, except from trying to find the author and tell him ? Is there a GitHub repository somewhere that would hold them ?

Plan B would be to upload a modified script, I guess...


texts: Changing date format of feature class with ArcPy?

texts:
I thought this was fairly straightforward; however I am receving NameError: name 'dt' is not defined.  What is being done incorrectly?  My orginal timestamp is m/d/yyyy h:mm:ss i would like to change it to %m/%d/%Y

 >>> rows = arcpy.UpdateCursor(aTable)
    ... for row in rows:
    ...     datetimeVal = row.getValue("last_edited_date")
    ...     formattedTime = dt.strftime(datetimeVal, "%m/%d/%Y")
    ...     print formattedTime
    ...     row.setValue("PickupDate", formattedTime)
    ...     rows.updateRow(row)
    ... del rows, row


texts:
I am trying to load S-57 file with pyQGIS with code below

self.canvas = QgsMapCanvas()
# self.canvas.useImageToRender(False)
self.canvas.show()
# Lay our widgets out in the main window using a
# vertical box layout
self.layout = QVBoxLayout(self.frame)
self.layout.addWidget(self.canvas)

file = QFileDialog.getOpenFileName(self, "Open S57file", ".", "S57files (*.000)")
fileInfo = QFileInfo(file[0])
# Add the layer
layer = QgsVectorLayer(file[0], fileInfo.fileName(), "ogr")
if not layer.isValid():
    print("Invalid Layer File")
    return

QgsProject.instance().addMapLayer(layer);
self.canvas.setExtent(layer.extent())
self.canvas.setLayers([layer])

But I am getting error

"ERROR 4: S57 Driver doesn't support update."
and file is canvas is blank although the layer is valid

I came across this question and also set the system path variable named s57_CSV referring to the required csvs folder but no difference.


texts: Performing multilayer intersection in QGIS?

texts:
This question on Multilayer intersection in QGIS was asked by "julsbreakdown" a year go. But I am a newbie and I am not really sure how it works with SQL as suggested in the responses. And I am using the same image posted by "julsbreakdown." Here is the link of the previous thread: Multiple layer intersection

Using QGIS (or its plug-in) or python script, how do I get the polygon (x) which is the intersection of four polygons?


texts: How to create point feature from longitude/latitude coordinates using ArcObjects and VB.net?

texts:
i Want to add a point feature to my feature class. an i want to enter the coordinates manually.

here is what i am trying to do

    pworkspaceedit = pFeatureWorkspace
    pworkspaceedit.StartEditing(True)
    pworkspaceedit.StartEditOperation()

    Dim point As IPoint = New Point
    Dim geometry As IGeometry

    point.X = 100
    point.Y = 100

    geometry = point

    pfeature = pFeatureClass.CreateFeature()

    pfeature.Shape = geometry

    MsgBox(pFeatureClass.FeatureCount(Nothing))
    pworkspaceedit.StopEditOperation()
    pworkspaceedit.StopEditing(True)

i have searched and the results are for copying from one feature class to another using icursor, creating lines from points, but havent found one for creating a point feature. all suggestions are welcome

there is no error in the code, a new empty feature is also created, but it has no attributes, MsgBox(pFeatureClass.FeatureCount(Nothing)) this i am using to display the no. of features after its created


texts:
I was wondering how Google Maps on my Android phone can locate me with a good accuracy when my GPS is disabled.

At first I though it was with the cell tower but I'm not sure it can achieve this precision and it is way better when the wireless is enabled.

If Wifi is ON and GPS is ON = Perfect location
If Wifi is ON GPS OFF = Still Perfect location
If Wifi is OFF and GPS is ON = location is not presisce.

Tested this with 5 -6 Phones

Any Idea why ?


texts:
I'm setting up an published arcreader file, and was searching for a way to add a date/time stamp to the printed map.

Is there a way to add a dynamic text box to a published pmf file that will show the date that the map was opened?


texts: QGIS version of Arc's "ungenerate" command?

texts:
I'm attempting to export the coordinates of the vertices of a vector layer in QGIS 2.8.9-Wien.  I need them to be in a specific format for a simulation program to load them.  The documentation for this simulation program says the following.

"The geometry file format is adopted from the “ungenerate” command of the ArcGIS geographic information system (other GIS packages can also use this format). “Ungenerate” exports the geometry of a polygon system as an ASCII text file."

The documentation then shows the format as so:

  1      2923.3222656      2084.7429199
  2917.0759277      2079.6772461
  2926.8261719      2089.6940918
  2929.4162598      2089.6560059
  2918.4470215      2078.3823242
  2917.0759277      2079.6772461
END
     3      2936.2717285      2087.4851074
  2929.4162598      2089.6560059
  2942.2133789      2087.4470215
  2942.0610352      2084.7048340
  2930.1779785      2087.5612793
  2929.4162598      2089.6560059
END
     4      2937.0334473      2081.5056152
  2942.0610352      2084.7048340
  2941.4516602      2076.0590820
  2934.9006348      2073.7739258
  2930.1779785      2087.5612793
  2942.0610352      2084.7048340
END
END

It is a ID number for each shape (1,2,3,...) then the X,Y coordinates of the vertices listed in pairs and finally "END" after all are listed.

Currently I can use the MMQGIS plugin to export a similar format which I reformat in excel.  I want to automate the process using the graphical modeler so I'm wondering a) is there a way to export from QGIS to get this format directly and b) is that method available in the graphical modeler.


texts: rgdal::writeOGR collision in column names results in NULL data

texts:
I already know (How to stop writeOGR from abbreviating Field Names when using “ESRI Shapefile” driver) that writeOGR is abbreviating names of columns longer than 10 chars. But there is one more thing. If some two column names are same in case-insensitive way, it will remove data from both/one of them!

library(rgeos)
library(rgdal)
library(sp)
line.data.df <- data.frame(a = 1, b = 2, Aa = 3, bB=4, aa = "hello", bb = "bug" , Depth.Stratum.max = 6, Depth.Stratum.min =3)
l <- rgeos::readWKT("LINESTRING(0 7,1 6,2 1,3 4,4 1,5 7,6 6,7 4,8 6,9 4)")
l.df <- SpatialLinesDataFrame(sl = l, data = line.data.df)
writeOGR(obj =l.df, dsn = ".",  layer = "test_abbreviation", driver="ESRI Shapefile", overwrite_layer = T)
l.shp <- readOGR(dsn = ".", layer = "test_abbreviation")
l.shp@data

With output

> library(rgeos)
rgeos version: 0.3-23, (SVN revision 546)
 GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246
 Linking to sp version: 1.2-5
 Polygon checking: TRUE

> library(rgdal)
Loading required package: sp
rgdal: version: 1.2-8, (SVN revision 663)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
 Path to GDAL shared files: /usr/share/gdal/2.1
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-5
> library(sp)
> line.data.df <- data.frame(a = 1, b = 2, Aa = 3, bB=4, aa = "hello", bb = "bug" , Depth.Stratum.max = 6, Depth.Stratum.min =3)
> l <- rgeos::readWKT("LINESTRING(0 7,1 6,2 1,3 4,4 1,5 7,6 6,7 4,8 6,9 4)")
> l.df <- SpatialLinesDataFrame(sl = l, data = line.data.df)
> writeOGR(obj =l.df, dsn = ".",  layer = "temp.samplings", driver="ESRI Shapefile", overwrite_layer = T)
Warning message:
In writeOGR(obj = l.df, dsn = ".", layer = "temp.samplings", driver = "ESRI Shapefile",  :
  Field names abbreviated for ESRI Shapefile driver
> l.shp <- readOGR(dsn = ".", layer = "temp.samplings")
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "temp.samplings"
with 1 features
It has 8 fields
> l.shp@data
  a b Aa bB aa_1 bb_1 Dpth_Strtm Dpth_Str_1
0 1 2  0  0 <NA> <NA>         NA         NA

You can see that strange and annoying abbreviation + case-insensitivity caused columns (especially my Depth.Stratum.min and Depth.Stratum.max) to have colliding names resulting in data loss.
Weird:

some letters in names are kept in UPPER case.
in case of columns (aa and aA) - data from one of them was preserved, in case of my columns Depth.Stratum.min and Depth.Stratum.max none contained original data

Is this behavior a bug or it's just something I have to live with?


texts: Detecting "short segments" in a File Geodatabase Feature Class

texts:
Using ArcMap 10.3 (without any extensions) is there a way to detect "short segments" in a feature class that is part of a file geodatabase automatically (without having to zoom to each feature to inspect manually)?

"short" is defined by a value the user enters. "segment" refers to the distance between successive vertices in a polygon or polyline.


texts: Why do I get different output from gdal2tiles vs. maptiler?

texts:
I am trying to generate map tiles using gdal2tiles for an area much larger than MapTiler can handle.  However, my tile caches are not showing up on the maps.  I checked gdal2tiles by running the same file through it and MapTiler - and the outputs come out differently.  Here's my gdal2tiles command:

  /usr/local/Cellar/gdal/1.9.0/bin/gdal2tiles.py -p mercator --s_srs
  3857 -z 1-9 -v  myFile.tiff

For instance, in the maptiler output, I get the file

  9/85/197.png

While gdal gives me different files, like

  9/85/313.png

The output should be identical, right? Is there something I'm missing?  My projection on the original file is definitely correct.  I have also tried gdal2tiles with 900913.  I am running gdal2tiles on MacOS X Mountain Lion, gdal 1.9.0.


texts: Clip polygon with nothing more than a point as starting point

texts:
Quite a straightforward question:
I have a point layer that consists of a few thousand points. These points are heterogeneously positioned in or very close to polygons that are water bodies. These polygons differ in size. I would like to cut the polygon based on the location of these points.
Please look at the example below. Blue polygon (water body) and black point is what I have right now. The purple, red or green lined polygon (or line) is a example of what the result can look like in order to clip the blue polygon. Which approach do you suggest?
Perhaps that a buffer around a point would help. The tricky part is that the buffer should match the two closest edges of the surrounding polygon. Or maybe to create a line to both sides?

Other suggestions?
Working with both QGIS and GIS (arcmap 10.0)


texts: How generate ID uniquely in two featureStores with GeoTools?

texts:
I have a PostGIS datastore in GeoTools as follow:

JDBCDataStore dataStore = ...

and I have 2 tables(or FeatureStore in GeoTools) as follow:

SimpleFeatureStore table1Store = (SimpleFeatureStore) dataStore.getFeatureSource("table1")
SimpleFeatureStore table2Store = (SimpleFeatureStore) dataStore.getFeatureSource("table2")

I can insert data in two tables. Inserted record ID generated automatically. I want generated ID in two tables be unique. For example if a record insert to table1 and id 22 generated, next insert to table2 generate ID 23.

My resolve is change default value of column id in table2 to nextval('table1_id_seq'), but when I want to insert in table2 occurs an error as follows:

org.postgresql.util.PSQLException: ERROR: currval of sequence "table1_id_seq" is not yet defined in this session.

What can I do to achieve this goal?


texts: Filter spatialite layer with 'not equal' in QGIS

texts:
I tried to filter my spatialite layer with the

"SPECIES" IS NOT NULL OR "SPECIES" != 'Corvus corax'

expression, but the result contains the lines where "FAJ" field is 'Corvus corax' (the NULL value withered correctly).

What's wrong with this expression? How can I filter all of the rows where "SPECIES" field is not null or not equal to Corvus corax?


texts: Why doesn't SQL sum(dollar) work in SQLITE?

texts:
I've got a sqlite table with real numbers, two decimal place precision, in a mydollar_column.  (There are no null values, but some of the values are negative numbers) For some reason, using sql on those values, always returns 0.00000.  I've tried this in Qspatialite (with QGIS 1.8) and Spatialite_GUI.

Here's an SQL example:   SELECT sum(mydollar_column) FROM myTable

If I use the sum function on an integer column, it works.  For example, I used sum(pkuid)on the primary key, and the sum() works as expected.

As a test, when I export the table to csv file, the values of my dollar_column are real numbers with two decimal place precision, for example "1234.56"  But Qspatialite and Spatialite GUI always shows the values as, for example "$1,234.56" I've also used the total(mydollar_column) function with the same unsuccessful results.

Does anyone know why I'm unable to calculate my dollar_column?
How does Qspatialite and Spatialite_GUI interpret the values of mydollar_column as dollars, when sqlite is typeless columns?