texts: Adding new layer to existing User Stored Display in ArcFM?

texts:
I am trying to help a coworker to add a brand new created layer in our GIS to the User Stored Display that is saved. The saved display existed previous to the new layer being added to the our GIS.

Is this as simple as 'adding data' while having the user stored display active and saving it again?

It would make sense to me that it is that simple but I can't play around with it since I started working later than the other person (who is old school) so I don't have the issue.

I have tried looking around and I can't find anything on it and I do was not able to just try adding the data because my coworker was working on something and they don't believe it is that simple. We are still using ArcMap 10.2. with ArcFM.


texts: Using training polygons to classify past imagery based on spectral reflectance in ArcMap

texts:
I am currently creating training polygons on ArcMap on 2020 imagery and classifying that imagery.  I want ultimately to compare it to historic imagery (I'll first do 2013 for example).  However, I want to be able to use the spectral data and not necessarily just the location of the training polygons on the historic imagery.
To be more specific, training polygons on 2020 imagery are used to classify what is bare dirt or what has vegetation, which might not be the same across all locations for both 2020 and 2013.  So is there a way I can use the spectral properties of training polygons to classify historic imagery.
I don't want the polygons to use the spectral reflectance of their location in 2013 but rather the spectral reflectance from 2020. I'm hoping this makes sense!


texts: MapServer vs QGIS Server for WebGIS with access control

texts:
I'm looking to set up a FOSS web GIS to allow sharing of large data files.  From online discussions I gather this can be done using a variety of software combinations, e.g.

PostGIS + MapServer + OpenLayers

or alternatively

PostGIS + QGIS Desktop + QGIS Server + QGIS Client.

Could people experienced in this area please detail any advantages and/or disadvantages that might be associated with these two options?

Are there any problems setting up either of these two systems on a 64-bit machine?

The end product webGIS would ideally give differential access via user name and password, with different selections of data available to different users.

Any advice appreciated.


texts: Multiple font markers in shapefiles in QGIS?

texts: I'm making a UK phase 1 vegetation survey map and require to show the dominant plant species in a given polygon.  In phase 1 these are shown as 2- or 3-letter codes e g Ri, Cys.  These codes need to be located very precisely and as far as I can see, the various labelling options in QGIS do not allow sufficient precision.  I am therefore trying to create point shapefiles for each dominant plant species.  Is there a way of using multiple letters in the "font marker" symbol layer type?


texts:
I need to create a map of my sampling site, and the journal requested to make an inset map showing where in the world this map is, this is my map code.

library(maps)
library(GISTools)

map('state', fill = FALSE, xlim = c(-125, -114), ylim = c(32.2, 42.5), xlab = "lon", ylab = "lat")
map.axes(cex.axis=0.8)

points(-121.6945, 39.36708, bg = "black", pch = 21)

maps::map.scale(x=-124, y=34, ratio=FALSE, relwidth=0.3)
north.arrow(xb=-116, yb=41, len=0.22, lab="N")

I have tried to find a way to make the inset map but no luck so far.


texts: Adding building/woodland elevation data to existing OS Terrain 5 DTM in ArcGIS Desktop

texts:
Would someone be able to give me the stepped process showing how to add building/woodland elevations to an existing DTM raster in ArcMap?

I seek this in order to prepare a surface model for a viewshed analysis.

I have run a viewshed analysis on a series of points using a DTM raster derived from UK OS Terrain 5. This creates the bare-earth model and is working well.

I now need to overlay this with buildings (at 8m) and woodlands (at 15m) from UK OS Open Map Local, and then add this elevation data to my existing Terrain 5 DTM. This is the bit I'm struggling with:

Adding the elevation data to the UK OS Open Map shapefiles for buildings and woodland. The method to do this would be to add an extra field in the attributes table and assign the 8m/15m elevation to this field. What field type/properties should the new field have?
Converting the shape file to a raster dataset using the Polygon to Raster tool on each raster created.
Merging with the existing Terrain 5 DTM using the Raster Calculator tool. How do you assign the cells with no trees/buildings a value of '0'?

I've seen the question answered for QGIS but not ArcGIS Desktop.


texts: How to export streets in a certain area from one shapefile to a new shapefile?

texts:
I have a shapefile with the roads from Europe and I want to export in a new shapefile only the roads from Greece.
I tried to delete the other roads but there are too many and I have no attribute to help as an example country.

the attribute is

  gid serial NOT NULL,
  sourceid character varying(50),
  exs numeric,
  notes character varying(254),
  roadid character varying(50),
  onme character varying(200),
  rtenme character varying(200),
  ntlclass character varying(200),
  fclass numeric,
  crgway numeric,
  numlanes numeric,
  lnewidthm numeric,
  rdwidthm numeric,
  axleloadmt numeric,
  totloadmt numeric,
  srftpe numeric,
  srfcond numeric,
  srfprep numeric,
  isseasonal numeric,
  curntprac numeric,
  gdwthrprac numeric,
  bdwthrprac numeric,
  speedlimit numeric,
  curntspeed numeric,
  gnralspeed numeric,
  isundrcstr numeric,
  cstwrketc date,
  graddeg numeric,
  sec numeric,
  hasshouldr numeric,
  hassidewlk numeric,
  drivside numeric,
  iselevated numeric,
  hasmedian numeric,
  opstatus numeric,
  length_km numeric,
  shape_leng numeric,
  geom geometry(MultiLineString,4326),

I use postgres and now start to use QGIS


texts: Problem with Sentinel-2 NDVI export to TIFF in Google Earth Engine

texts:
I'm trying to calculate and export NDVI rasters extracted from Sentinel-2 using the script:

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

Source code: Export images in google earth engine with date name

When I run it for coordinates

  [27.493905, 53.73073],
  [27.493905, 53.990091],
  [27.849634, 53.990091],
  [27.849634, 53.73073]

it works properly and I get correct NDVI values (-1 - 1).

But when I run this script on  coordinates

  [47.11727, 54.71504],
  [47.11727, 54.942041],
  [47.411848, 54.942041],
  [47.411848, 54.71504]

it still works and show the first calculated NDVI in GEE but after I export rasters on local computer from Google Drive I get strange values of NDVI and can't visualise NDVI in GIS...

Low : -3,40282e+38
High : 3,40282e+38

How can I solve this issue?


texts:
I seek advice for calculating presence/absence of features in QGIS.  I have a polygon dataset of boundary outlines and a number of other datasets which I want to compare it against for presence/absence of features within the main boundary dataset.

Can I do this using the vector datasets as is or will I need to convert then to a grid/raster beforehand?

My example is I have OS administrative boundaries an OS open green spaces data, I would like to know whether an area of open green space is present or absent from each boundary polygon. I was thinking about using the Join attributes by location and intersect, but don't know if that is flexible enough to give a simple yes/no answer for each boundary polygon.


texts: Dragging and dropping from Catalog Favorites of ArcGIS Pro

texts:
I'm slowly transitioning from ArcMap to ArcGIS Pro (2.3.0) and I'm facing a quite frustrating UX problem:

I can't drag and drop layers from the Catalog pane to the map if it comes from the Favorites menu. I have to right-click on it > Add to current map.

Oddly, drag and drop works with layers from folder connections created in the Project tab but those folder connections are linked to the project so you have to recreate them for every project.

It might not seem like a big issue but it's one of the little things that bother you when you work all day with this software. I just wanted to check if this is the same for any of you, or if someone found a way to solve it (maybe with a newer version?).


texts: Using ArcPy and Python to export Table View as PDF format?

texts:
I am using ArcPy to select from a layer and display the result on map as a table view.

Now I want to export the table view as a PDF format.

Anybody know about it?

If it is beyond ArcPy to export as PDF directly, with intermediate format also welcome.


texts:
I'm using Leaflet for the maps of my website. So far I was using maps of Mapbox (they're cool), but they obviously charge for them for website with 100,000+ maps per month.

I'm trying to find a provider of tiles which no charge for their tiles. Initially, I do it with OpenCycleMap, and substituted

https://{s}.tiles.mapbox.com/v3/{id}/{z}/{x}/{y}.png

with

http://{s}.tile.opencyclemap.org/cycle/{z}/{x}/{y}.png

Two questions:

1) Does OpenCycleMap provide free tiles with no restriction of maps per month?

2) Does anybody know any other free provider?


texts: Getting number of pixels from raster that fit a certain criteria

texts:
I have downloaded many tiles of MODIS' MCD64A1 Burnt Area GeoTIFF data per year and am trying to count the number of pixels that are not -2, -1 or 0 (to calculate approximate total area burnt each year).
I have tried direct approaches to make statistical analysis, such as using geopandas and pandas on Python, but to no avail, since it returns errors that say I must allocate 90 plus pebibytes of memory (probably due to the data being nation-wide).
What are some ways to extract the number of pixels from multiple rasters that are not certain values?


texts: How to determine if one polygon is contained in another?

texts:
Using the Net Topology Suite, what is the easiest way to determine if one polygon is entirely contained in another?

Given two independent polygons, A and B, they may overlap, or they may not. If they do, polygon A may be entirely contained in polygon B or polygon B may be entirely contained in polygon A.

I'm looking for an answer to the simplest case, is polygon A entirely contained within Polygon B.

Is the simplest test simply to see if each point of A lies within B?


texts: Calculating length of line using British National Grid projection?

texts:
I need to calculate the length of a line in (KM) in QGIS using the British National Grid as the Projection.

Can somebody explain to me how to do this?


texts: Relating a table to a point layers and exporting results in model builder?

texts:
Here are the data that I'm working with:

Table:

Standalone table from an OLE DB connection with a MS Access db
10,800 records with non-distinct account IDs
Divided into 300 different categories

Layer:

Spatial point layer
1,700,000 records with distinct account IDs

I ultimately need to create a compiled feature class made up of points representing each of the 10,800 account ID records in my Access database. The idea is that some of the points will be duplicated and redundantly stacked in this final layer — I could do one giant relate, but I need the many-to-one relationship symbolized.

The Account IDs are only duplicated between the 300 categories, i.e., Category1 has distinct account IDs, but you could see some of the same IDs from Category1 in Category2.

Here is my manual workflow:

Select just the records for Category1 in the Table
Select related points from the point layer and create a point layer selection
Join the data from the connected table with the point selection and export as a new feature class
Merge/Append with ultimate compiled feature class

Doing this for all 300 categories seems tedious. I've experimented with "Make Table View" and "Table to Excel" but it doesn't look like I can relate or join anything in the ModelBuilder and have it spit out a new selection feature class. Plus, it would be great if I could have it iterate through the relates based on all of the 300 unique category values.


texts: Unable to load PostGIS vector layer into QGIS using python cookbook recipe

texts:
I'm having trouble loading a PostGIS layer using Dufour (from osgeo4w), using the guide at http://documentation.qgis.org/2.0/en/docs/pyqgis_developer_cookbook/loadlayer.html. I just get the 'invalid layer' message. My parameters appear to be perfect.

from PyQt4 import QtGui
from PyQt4 import QtCore
from qgis.core import *
from qgis.gui import *
from ui_townplanning import Ui_TownPlanning
# create the dialog for zoom to point

class TownPlanningDialog(QtGui.QDialog):
    def __init__(self, iface):
        QtGui.QDialog.__init__(self)
        # Set up the user interface from Designer.
        self.ui = Ui_TownPlanning()
        self.ui.setupUi(self)
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
        self.addPostGISLayer('myhost', '5432', 'glis', 'user', 'pass',
                             'tpadmin', 'planning_application_boundary_validation',
                             'wkb_geometry')

    def addPostGISLayer(self, host, port, dbname, username, password,
                        schema, table, geom_col):
        '''Adds a PostGIS table to the map'''
        uri = QgsDataSourceURI()
        uri.setConnection(str(host), str(port), str(dbname), str(username), str(password))
        uri.setDataSource(str(schema), str(table), str(geom_col), '', str(key_col))
        vlayer = QgsVectorLayer(uri.uri(), str('test'), 'postgres')
        # QgsMapLayerRegistry.instance().addMapLayer(vlayer)
        if not vlayer.isValid():
            msgbox('Layer not loaded', uri.host(),
                   uri.database(), uri.port(), uri.username(), uri.password(),
                   uri.schema(), uri.table(), uri.geometryColumn())

def msgbox(*args):
    a = [str(i) for i in args]
    text = '\n'.join(a)
    qmsg = QtGui.QMessageBox()
    qmsg.setText(str(text))
    qmsg.exec_()

Code now works after correction


texts: How to fix problems opening JPEG2000 in QGIS Mac version 1.8?

texts:
I am running QGIS v. 1.8 on both my Windows and Mac partitions; Every time I try to open a jpeg2000 (jp2) raster on the mac version the program locks up and I need to force quit.

Does anyone have any ideas how to fix this problem?


texts: Using expressions in generic QGIS Print Layout inserted tables

texts: I would like to use expressions in tables added in QGIS Print Layout using the Add Fixed Table option, as shown in the screenshot attached the expressions are instead displayed as regular text in the table's cells


texts: ArcMap - Making DEM using GPS vector data

texts:
I am new to GIS and GPS. I am going to make a digital elevation model of my school's campus using a trimble gps unit. I know the accuracy is not the best but it is just for a school project.

My question is how exactly should I go about doing this? I went out in the field today and acquired a bunch of position points which have elevation (I set the gps to log my position every second). I made that trimble data file into a .shp file and put that into arcmap. I can view all of the points and the elevation in the attribute table.

How to I make a digital elevation model from this data? I tried doing raster interpolation like kriging and IDW but I am not really sure what I am doing. What is a TIN? I looked TIN up and apparently it uses Vector data? Shouldn't I just use my vector data?


texts: Georeferencing raster using world file in ArcMap?

texts:
I have downloaded a raster (jpg) and associated world file (jpgw) from https://qimagery.information.qld.gov.au

I have turned on the option in ArcMap to use the world file to define the coordinates of the raster.

When I bring the file into ArcMap I get the warning that there is an unknown spatial reference.

I have tried editing the spatial reference for the raster in ArcCatalog for the file for a number of common protections used in Australia. I always get the same error. Perhaps I'm just tried the wrong coordinate systems?

Interestingly, when I bring in a raster and world file sourced from nearmap instead, it locates the image in the correct place, without me having the set a spatial reference for the raster, but still gives the "unknown spatial reference" warning.

I've attached the image and world file I'm having trouble with under the following links:

https://drive.google.com/drive/folders/1GjeBcFDQmH7coYJtq7CQlZdB9TgzJmeI?usp=sharing


texts: Converting TIFF (DIMAP) to .xyz file format with GDAL

texts:
We are attempting to convert a .tif file into .xyz file.

We have been using gdal_translate on a linux box and this is the line we've used:

gdal_translate -strict -of XYZ ABC.tif ABC.xyz -mo ABC.dim

The file has three bands: NIR, red and green respectively and we are trying to extract these values.

The dim file is just the meta data details describing the .tiff file.

I currently get a segmentation fault when I run this line, can anyone tell me how to fix this?


texts: Appending Attribute Table of multiple feature class with same points/ pointID

texts:
I used extract by points function to extract values from multiple rasters using a point feature class. The output is a couple hundred point feature classes with values associated with those rasters.

I would like to append all the attributes table of feature class (they all have same POINTID) so that I have all the information I want in a single table.

How do I do this?

See attached screenshot:

I would like to append RASTERVALU of each feature class (SRAD_KRIGXXXXX) to "FOOTPRINT" attribute table (maintaining its file naming structure).

Such that the FootPrint feature class will have field names SRAD_KRIG2007001, SRAD_KRIG2007002, SRAD_KRIG2007003, etc. with values from their RASTERVALU field.

I was wondering if python scripting is what I need.


texts:
This may be a simple question but can't seem to find an answer. I've got a raster with 5 attributes, and a SpatialPointsDataFrame. I would like to count the number of points that fall within each attribute of the raster (i.e. how many points fall within attribute 1, how many within 2, etc...).

My raster looks like this:

 > ratify(r.class)
 class      : RasterLayer
 dimensions : 168, 432, 72576  (nrow, ncol, ncell)
 resolution : 0.8333333, 0.8333333  (x, y)
 extent     : -180, 180, -56, 84  (xmin, xmax, ymin, ymax)
 crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
 source     : memory
 names      : layer
 values     : 1, 5  (min, max)
 attributes :
 ID
 1
 2
 3
 4
 5

How can I do this?


texts: What are the differences in raster formats between Civil3D, QGIS, and ArcGIS?

texts:
Civil3D uses a .dem / .tif format which calls many files from other sources invalid. Esri products seem to have richer metadata for TIFF rasters and will add data to TIFF files on disk. QGIS uses GDAL but doesn't automatically add any metadata.

I have found how to add pyramids, for example, in QGIS, but when I ingest a raster TIFF that I have exported from CloudCompare, for example, the metadata is basically non-existent. I would like to be able to route CloudCompare raster exports into C3D ultimately.

What are some of the specifications that determine the way GeoTIFF is implemented in these different platforms?


texts:
The county I work for currently has an ArcGIS online web map on their public works website that shows our capital improvement projects as polygons with a PDF linked in the popup showing project info and updates: http://www.yorkcounty.gov/380/Public-Works
Right now every time there is an update, we find the project in an Excel spreadsheet, save the sheet as a PDF, and swap out the old PDF attachment for the new one in the ArcGIS online layer. This has become labor intensive so we are looking for a way to streamline this processes into as few steps as possible.

What we are looking for is a way to automatically add updates to a project document and have them linked to their respective projects in the map?


texts: QGIS- Can't find fitting CRS (rotated regular grid)

texts:
I am using QGIS 2.18.15. I got raster layer data’s for wind speeds in Europe.
 (ftp://ftp-cdc.dwd.de/pub/CDC/grids_europe/monthly/wind/MIKLIP_DECREG/v002/)
The data description describes the projection as an:

  "rotated regular grid, virtual North Pole at 39.25 N, 162.00 W (rotated coordinates)"

I can add the raster to QGIS but I do not find the fitting CRS. I looked for a possible Coordinate Reference System online but could not find one. Do you have an idea which CRS might suit in this case?


texts: display openlayers lat/long map bound on page

texts:
I want to display the bounds of my current extent on the page. I am not too great with javascript. I know it must be related to this:
map.getLonLatFromPixel(new OpenLayers.Pixel(bounds.left, bounds.bottom))

Can't seem to get it. Any help would be great! Thanks.


texts: Creating Loop to replace text in page layout in multiple MXDs instead of one MXD at a time?

texts:
I found some code on a website that will perform a search and replace on page layout text elements. Currently you can only do it one MXD at a time. I was hoping to change the code so I can do a search and replace on a folder with a bunch of MXDs in it.

The script is below. I know I need  a loop, just not sure how to do it; never programmed before.

import arcpy, string, os

#Read input parameters from script tool
Path = arcpy.GetParameterAsText(0)
oldText = arcpy.GetParameterAsText(1)
newText = arcpy.GetParameterAsText(2)
case = arcpy.GetParameter(3)
exact = arcpy.GetParameter(4)
outputMXD = arcpy.GetParameterAsText(5)

try:
    #Referent the map document
    mxd = arcpy.mapping.MapDocument(mxdPath)

    #Find all page layout text elements
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if exact:
            if case:
                if oldText == elm.text:
                    elmText = elm.text.replace(oldText, newText)
                    elm.text = elmText
            else:
                if oldText.upper() == elm.text.upper():
                    elmText = elm.text.upper().replace(oldText, newText)
                    elm.text = elmText
        else:
            if case:
                if oldText in elm.text:
                    elmText = elm.text.replace(oldText, newText)
                    elm.text = elmText
            else:
                if oldText.upper() in elm.text.upper():
                    elmText = elm.text.upper().replace(oldText, newText)
                    elm.text = elmText
    mxd.saveACopy(outputMXD)

    del mxd

except Exception, e:
    import traceback
    map(arcpy.AddError, traceback.format_exc().split("\n"))
    arcpy.AddError(str(e))


texts: Add features to vector layer from WFS protocol read response

texts:
I'm trying to add the features I get from a WFS query using the WFS protocol. But something doesn't work.

I have looked at other posts with people having similar problems, but their solution doesn't fit with my problem.

I have tried:

var gmlReader = new OpenLayers.Format.GML({ extractAttributes: true });
var features = gmlReader.read(evtResponse.responseText);
markerVectorLayer.addFeatures([features])

But then I receive a error:

TypeError: a is undefined

...b.push(d)}return b},parseFeature:function(a){for(var b="MultiPolygon Polygon Mul...

OpenLayers.2.13.1.js (line 380)`

When I use console.log to check the features it looks like a perfect FeatureCollection object.

If I loop the response and create polygons from the coordinate points it works, but the performance isn't functional.

What am I doing wrong?


texts:
I need to make a query for my local OSM server which returns a table with 2 columns: user.display_name and the total of nodes of the ways with tag k='waterway' created, edited or deleted by that user.

I think it should be something like:

SELECT
    users.display_name,
    COUNT(current_way_nodes)
FROM
    users
INNER JOIN changesets ON users.id = changesets.user_id
INNER JOIN current_ways ON changesets.id = current_ways.changeset_id
INNER JOIN current_way_tags ON current_ways.id = current_way_tags.way_id
INNER JOIN current_way_nodes ON current_ways.id = current_way_nodes.way_id
WHERE users.display_name = 'foo' AND current_way_tags.k = 'waterway';

But this query gives me the following error:

ERROR:  column "users.display_name" must appear in the GROUP BY clause or be used in an aggregate function
LINE 2:     users.display_name,

If someone need to have a look at the schema, it's available here: https://wiki.openstreetmap.org/w/images/5/58/OSM_DB_Schema_2016-12-13.svg

I appreciate any help.


texts: Get Spatial Table Information Oracle Spatial

texts:
How can I get the tablename, Srid and GeometryType of all spatial tables in Oracle(I use 11g)?

In postgis, I can do something like

select f_table_name,srid,type from public.geometry_columns order by f_table_schema, f_table_name;

Does anything similar exist in Oracle?


texts: Changing coordinate system with rasterio affects dataset opening

texts:
I have troubles changing coordinate system of a raster with rasterio. Here is my code:

Here I create my .tif (it's a RGB image)

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

    with rasterio.open(result_path + "\\" + "test_infra_{}.tif".format(product_group), 'w', **meta) as dst:
        dst.write(bandstack_infra)
        dst.close()

Then I re-open it and change the meta and create a new .tif with the new coordinate system, according to the piece of code given in the docs (https://rasterio.readthedocs.io/en/latest/topics/reproject.html).

    dst_crs = 'EPSG:4326'
    with rasterio.open(result_path + "\\" + "test_infra_{}.tif".format(product_group)) as src:
        transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height})

    with rasterio.open(result_path + "\\" + "test_infra_{}_wgs.tif".format(product_group), 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

Then I get the error:
rasterio.errors.RasterioIOError: Dataset is closed: --> path to the original data, open as src.

I understand that the new dataset cannot be created if the original one (open as src) is not open. However, I didn't close it (with src.close())!
Why is it closed and how can I solve this problem?


texts: How can I classify GPS points with road type?

texts:
I have a road database in .osm.pbf format; I want, for a long list of GPS points, to classify whether they fall on (with a margin of, say, 20m)

Highway: motorway, trunk, motorway_link, trunk_link
Intercity: primary, secondary, tertiary, primary_link, secondary_link, tertiary_link
Other

I'm open to use any tool that will do the job - arcpy, spatialite, postgis, others - but I'm not deeply familiar with any of these.


texts: ST_Difference in PostGIS for cutting out areas of polygons

texts:
I am trying to understand SQL commands in PostGIS.
I am trying to perform an operation which would be completed with the use of the Difference Operator in QGIS but I have large numbers of polygons so I need to manage this in PostGIS instead.
Essentially, I have some circular polygons. Let's call those buffers and I want to cut out some areas from those buffers. Let's call the areas to be cut out "blanks". The output should be the buffers minus any areas where they overlap with the blanks.
I've looked around for some suitable code to re-use and have found the code below.  This works, after a fashion, BUT it doesn't give me the buffers which do not intersect with any blanks, which I also need. The problem is, I don't understand how this is working.
Why does JOIN use the expression FROM buffers b JOIN blanks a i.e. why isn't it FROM buffers JOIN blanks?
Why do we need an ST_INTERSECTS as well an ST_DIFFERENCE?
Shouldn't ST_DIFFERENCE on its own be enough?
Why is it so complicated to create what is presumably a regularly used operation?
I am new to SQL and these may be obvious questions. I've worked my way through as many tutorials as I can but they all seem to be quite basic or too advanced for me to follow. The second things progress beyond simple expressions I get lost.
CREATE TABLE buffers_minus_blanks AS
WITH temp AS
(
  SELECT   b.id, st_union(a.geom) AS geom
  FROM     buffers b JOIN blanks a ON st_intersects(a.geom, b.geom)
  GROUP BY b.id
)
SELECT b.buffer_id, st_difference(b.geom,t.geom) AS newgeom
FROM buffers b LEFT JOIN temp t ON b.id = t.id;


texts:
I'm trying to write a timelapse function in Google Earth Engine, but the standard JavaScript setTimeout (which I'm using to set a delay between frames) does not appear to work, as it complains about it not being defined in this scope (I imagine the API does not access the DOM the same way a webpage would). I have a video export function that works, but I'm also trying to render something like this in-engine.

Is there some way to essentially tell GEE to wait, or something analogous to the setTimeout js function that I'm overlooking?


texts: How to add a button to the Qdockwidget Layers in Qgis 3 with PyQgis

texts: I want to add a plugin for styles to the layers panel. So I am trying to find a way to add a button to the toolbar of the layers-panel. Is this possible? If not is it possible to add a menu-item to the right-click menu of that same panel.


texts: What free gis is on android that leads in the area? -- just name top 1-2

texts:
similar kind of q of to What are some Free and Open Source GIS Desktop packages

include evidence such as 4.8-5 rating out of 5, etc.


texts: PostGIS views vs WFS Filtering operations on PostGIS table

texts:
In what use cases should I prefer writing PostGIS views for filtering datasets versus filtering with WFS filter encoding operators and vice versa?  In both cases I would be serving the data via GeoServer through the PostGIS table or view.

The main pro I can think of for the WFS filter approach is that the query can be encoded easily through the front end according to the filter specifications without needing to add extra logic to the backend code to handle the PostGIS queries.  Also, how does the WFS filter fare performance wise against PostGIS?


texts: QGIS 2.14 Couldn't load PyQGIS. No module named Qsci

texts:
When I open QGIS Version 2.14 Essen, I get the following error message:

Couldn't load PyQGIS.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/PROGRA~2/QGISES~1/apps/qgis/./python\qgis\gui\__init__.py", line 26, in
    from qgis._gui import *
ImportError: No module named Qsci
Python version:
2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

QGIS version:
2.14.0-Essen 'Essen', 98e9610

Python path:
['C:/PROGRA~2/QGISES~1/apps/qgis/./python.
...
...
...


texts: Joining a shapefile attribute table to a CSV file pyqgis3.0

texts:
I would like to join a shapefile to a CSV file in PyQGIS 3.0. So, when I run the following code, it executes without errors but the join action does not succeed. What is wrong with the following code?

layer=iface.activeLayer()
lien_csv="D:/Ghaleb/Projet_SIG/PFE SIG/FICHIERS_SHP/Essai.csv"
infoLyr =QgsVectorLayer(lien_csv,'Table_CSV' , "ogr")
QgsProject.instance().addMapLayers([infoLyr])
infoLyr.isValid()`
info = QgsVectorLayerJoinInfo()
info.setJoinFieldName = "CLE_CSV"
info.setTargetFieldName = "CLE"
info.setJoinLayerId = infoLyr.id()
info.setUsingMemoryCache = True
info.setJoinLayer(infoLyr)
layer.addJoin(info)
lien_save="D:/Ghaleb/Projet_SIG/PFE SIG/FICHIERS_SHP/Vanne_Joint.shp"`
QgsVectorFileWriter.writeAsVectorFormat(layer,lien_save, "iso-8859-15",QgsCoordinateReferenceSystem(32632, QgsCoordinateReferenceSystem.EpsgCrsId), "ESRI Shapefile")
joinedLyr = QgsVectorLayer(lien_save, 'Joined' ,"ogr")
QgsProject.instance().addMapLayers([joinedLyr])


texts: Making question conditionally required in Survey123 form?

texts:
When building a survey in Survey123, you can choose whether or not to make a question required.  But how do you make a question conditionally required?  For example, you have the following questions:

Do you like GIS?  (Yes or No?)
Why not? (Required if the answer to Question # 1 is "No")


texts: I have downloaded a CSV data set from Uber, but it gives locations in terms of Open Street Map IDs (and start and end node IDs) which Carto is not recognizing when I import. Is there some way to get Carto to recognize the locations through reorganizing my data?


texts: Exporting Multiple Variable Inputs to CAD in ModelBuilder?

texts: I am trying to create a tool that clips (X) layers from an area, and creates a CAD file with the layers. I'm using inline substitutions to fill in the names for the layers, (X), which are then supposed to be exported as CAD files. When I run the tool the way it is shown, all of the clips are created and sent to the appropriate folder, but only one CAD file is created. I'm assuming the issue is within the naming of the CAD Output, where a CAD file is created with the %Input Features% substitution, then each different layer that is exported is just overwriting the last due to the same substitution name.


texts: Creating Point Grid over selected polygon using OpenLayers 3?

texts:
I am trying to create an equally spaced grid of points. But these points must only be on land.

I was using the polygon overlay of countries to see if I can use that as a mask.

I have also tried using graticule to see if its possible to point points in a grid from this.

Documentation is very lacking on this. I can see how to do this using pointGrid() in OL2 but obviously being outdated this isn't going to be practical.

Can anyone give any pointers on how to over come this or where to start importing the pointGrid() function to solve the problem.


texts: How to make a displayed geoJson property a link to another property in a Leaflet popup

texts:
I have markers on my map created by geoJson data that represent restaurants, bars, etc.... Each marker has 2 properties; 1)Name - name of the site, bar restaurant, etc... 2)Link-a link to the establishment's website. Currently, I have added a popup that displays the Name. My goal is to make the name based on  into a hyperlink based on the link property that the user can click on and open corresponding the website. How could I best achieve this?

var layerGroup = L.geoJSON(data, {
  onEachFeature: function (feature, layer) {
    layer.bindPopup(feature.properties.Name +'
'+feature.properties.Link);
  }
}).addTo(map);

My best guess would be something like this: layer.bindPopup("" + feature.properties.Name + "")


texts: Merge shapefiles using ogr2ogr and retrieve the attributes of both

texts:
I have two shapefiles that I want merge (overlay/intersect) having in the destination the attributes of the two sources file.

ogr2ogr file_merged.shp tl_2010_us_uac10_conus_albersUA.shp
ogr2ogr -update -append file_merged.shp tabblock2010_10_pophu.shp -nln file_merged

In this way I obtain the polygons of the two shapes, but I lose the attributes of the tabblock2010_10_pophu.shp.

My final aim is then make an SQL statement where I can retrieve polygon in common in the both shapes.

Any suggestions??


texts: How to memorize temporary layer from processing.runalg() in pyqgis?

texts:
I created some geoalgorithms with QGis and I need to use them in a Python script. By saving the layer obtained in a file having a specific path, I do not have
problems.
For example:

processing.runalg("modeler:esplodilinee", "buildings", "C:\prova\sides.shp")
lines = iface.addVectorLayer("C:\prova\sides.shp", "poligon sides", "ogr")

But I would not use a file, but rather a temporary variable, in order to re-use it within the specific script.
The description of the algorithm suggests to me that:

ALGORITHM: esplodilinee
        buildings
        OUTPUT_ALGQGISEXPLODELINES_1

So, I thought to do something like this:

OUTPUT_ALGQGISEXPLODELINES_1 = processing.runalg ("modeler: esplodilinee"
"buildings", None)
lines = processing.getObject (OUTPUT_ALGQGISEXPLODELINES_1)

But I get a mistake and not the desired results.
Someone can help me?


texts:
So I am at a bit of a loss here, trying to get the distance in meters between a point stored in my table and a set of coordinates

locations

id    |   name  |    location(point)
-------------------------------
1        church     (52.632466,4.743884)

I'm using the following query for this

SELECT ST_Distance(
        ST_Transform(ST_GeomFromText((select location from locations where id = 1) , 4326),2163),
        ST_Transform(ST_GeomFromText('LINESTRING(-72.1260 42.45, -72.123 42.1546)', 4326),2163)
    );

And all I get is an error that the function does not exist, I have tested the function from the example on the documentation and that does work.

ERROR:  function st_geomfromtext(point, integer) does not exist
LINE 2:    ST_Transform(ST_GeomFromText((select location from locati...
                    ^
HINT:  No function matches the given name and argument types. You might need
to add explicit type casts.
SQL state: 42883
Character: 37

Any help in the right direction would be great.


texts: Understanding ArcGIS Geostatistical kriging?

texts:
We use soil samples data for interpolating soil fertility parameters using ArcGIS geostatistical Kriging. While doing kriging, the resultant map values ranges are not matching with our input samples range. Suppose available phosphorus of analysed soil samples data vary from 30 to 75 ppm (for 500 samples). After prediction, map shows range between 35 to 60ppm. Therefore we can't classify upper concentration, above 60 - 75 ppm.

Is it because of prediction modelling by calculating statistical relationship of input points?

What is the actual principle behind this?


texts:
I am currently taking my first steps in playing around with programming a plugin in python 2.7 for QGIS 1.8. While my plugin does what it should, I would like to make it more userfriendly. At the moment, the user is expected to choose an input- and an output-file, which are then processed. In case the user forgets to choose one (or both) of the files, the plugin is displaying an information window - and then quits completely after the user hits "OK".

if src_filepath =="" and dst_filepath =="":
    QMessageBox.information(None, "DEBUG:","Please specify files!")`

However, I would rather like to reload the main plugin-window with the "choose file"-dialog. I suspect that the solution is rather easy, but I was not able to invoke this window, so any help is appreciated.
The interface code derives more or less directly from the "Plugin Builder" in QGIS.


texts: Can I change the pixel sizes in a world file?

texts:
I start with a dataset that has projection WGS84 and pixel size 0.000274658203125 degrees/pixel. I run gdal_translate to convert it to a PNG:

gdal_translate dataset.tif output.png -of PNG -ot byte -scale -strict -co worldfile=yes

Running gdalinfo output.png shows that the pixel size has lost precision and is now 0.000274658200000.

However, I can edit the world file (output.wld) and change the X and Y pixel sizes to 0.000274658203125 and -0.000274658203125, respectively. Now, running gdalinfo output.png shows the proper pixel size.

So, my questions are:

Is pixel size only controlled by the world file and can therefore be changed simply by editing the world file?
If so, would it be better to construct my own world file in the first place?


texts:
In order to calculate TWI, I had to generate a catchment area raster, but I don't understand what the values of pixel represent. Any clarification?
The DEM is projeted in EPSG : 2154, so the unit is meter
I used the SAGA 'catchment area' tool


texts: How could I rate roads for scenic drives?

texts:
relative GIS newbie here.

I've been working on a personal project to rate roads on how nice a ride they'd be on a motorcycle.  I wrote a Python program to work through road data in the TIGER dataset (PostGIS) and apply a small number of metrics:

Road variance (twistyness)
Road elevation variance
Does the road pass by a body of water
Does the road go through a park

The results of applying the metrics have been a good start, but they are far from perfect.  So I'm trying to decide on some new metrics to apply.

Check land coverage around road.  Riding along rolling hillsides is better than a forest which is better than an industrial park.
Road condition (paved or unpaved)
Number of stop signs, street lights encountered
Utilize viewshed analysis to check if a body of water is actually visible from road
Identify mountain views and use viewshed analysis to see if mountain can be seen from road (think Great Smokey Mountains)
Apply historical traffic data

I need the advice of some GIS professionals.  Do these sound feasible, or even make sense?  Can you think of any other things I should try?

Most importantly, where can I get data for the valid ideas?


texts: How can I allocate a $rownum in reverse order in QGIS?

texts:
How can I allocate the $rownum value in reverse order in QGIS?

I have a created a grid layer to use in an  atlas and the ID allocated is from top left to bottom right - which is what I want if I want to create a $rownum attribute to use as the map number in my atlas (e.g. 1 - 78)

But I have selected only those grid polygons that have relevant data beneath them for my maps, so I have created a grid layer that is a subset of the original. In this subset grid, the ID numbers have been reversed - that is, they go from bottom right to top left so if I create the $rownum field the numbers will be in that reverse [new] order. (e.g. 45 - 1).

So I want to allocate that $rownum attribute in reverse order so that the atlas maps read from top left to bottom right.

Can anyone help?


texts: Error when running Extract by Mask Script

texts:
I am using ArcGIS 10.3 Desktop to write a script tool to, among other things, extract raster data using a mask.

I keep getting these errors:

  ExecuteError: ERROR 000875: Output raster: C:\Users02183940\Documents\Thesis_Data\soil_data\Calculated_Rasters\IL_rasters.gdb\Extract_ph_s1's workspace is an invalid output workspace.
  ERROR 000581: Invalid parameters.
  Failed to execute (ExtractByMask).

My code:

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

ebm_in_raster = r'C:\Users\a02183940\Documents\Thesis_Data\soil_data\Calculated_Rasters\ph_statsgo'
ebm_mask_data = r'C:\Users\a02183940\Documents\Thesis_Data\Cropland_Data_Layer\noda_1_10-15'

arcpy.CheckOutExtension("Spatial")
env.workspace ='C:\Users\a02183940\Documents\Thesis_Data\soil_data\Calculated_Rasters\IL_rasters.gdb'

outExtractByMask = ExtractByMask(ebm_in_raster, ebm_mask_data)
outExtractByMask.save('C:\Users\a02183940\Documents\Thesis_Data\soil_data\Calculated_Rasters\ph_in_crl')

What are the requirements for a valid output workspace for this tool? I don't think it is an issue of raster name length or use of unsupported characters as I've tried a few different things and none of them have worked.


texts: Create serial field following alphabetical order with QGIS?

texts:
I have a field with names, a string field, and I want to create another one with serial content following the alphabetical order.

I have this:

--------
|Field1|
--------
|  B   |
--------
|  C   |
--------
|  A   |
--------
|  D   |
--------

And I want to reach this:

-----------------
|Field1| Field2 |
-----------------
|  A   |    1   |
-----------------
|  B   |    2   |
-----------------
|  C   |    3   |
-----------------
|  D   |    4   |
-----------------

How can I reach it using the fiel calculator of QGIS?

Any idea???


texts:
I work with 64-bit Bentley WaterGEMS and it uses 64-bit Microsoft Access Database Engine. Opening XLSX files in ArcMap (via Add XY Data ...) needs 32-bit Microsoft Access Database Engine. Both 32 and 64-bit Microsoft Access Database Engine can not be installed simultaneously.

How can I do with this problem? Does ArcMap have 64-bit version?


texts: Populating ArcMap column with inputs from other colum, shortening original text?

texts:
I need to keep the first part of all inputs and get rid of the last part, so that i can later make a join based on the new attribute.

Edit: found solution by using simple python script !MY_FIELD![:-9], which put in field calculator removed the last 9 characters.


texts: Join a set of sequential points with individual line segments QGIS

texts:
My data set contains radar plots for flights on approach to an airport. For each unique flight I have hundreds of points as they are recorded every four seconds. Using Points2One I can join the points to form individual flight lines for each unique flight. In using Points2One by its nature it only creates a single line and a single record entry in the attribute table hence I'm loosing the data granularity of altitude and speed at each four second segment of the flight. Is there a way that I can create a single flight path that is composed of individual line segments, where each individual line segement would have the same record identifier (Operation No) and attrbutes as either the starting or ending point node of that line segement.

Example radar plot data

So with the data sample I'd like to end up with eight line segments that all have the Operation No. 16173263 and each segment would have the time delay, Latitude, Lon, Altitude & Velocity attributes of the starting point.

The radar plots for a single day contain approximately 1200 different flights and over 100,000 data points.

I have no python or coding experience and so far all QGIS experience consists purely of using menu driven functions.

Operation No, Time Delay (seconds), Latitude, Longitude, Altitude (ft), Velocity (knots)
16173263,0,42.005191,-81.987368,12444,0
16173263,1,42.003938,-81.986167,12418,328.5
16173263,5,41.999006,-81.981324,12303,326.6
16173263,9,41.994118,-81.976394,12182,326.6
16173263,13,41.989301,-81.971414,12054,324.6
16173263,17,41.98452,-81.966409,11919,322.7
16173263,21,41.979784,-81.961404,11778,320.7
16173263,25,41.975057,-81.956437,11637,320.7
16173263,29,41.970348,-81.95152,11490,318.8
16173264,0,41.512822,-82.082972,5331,0
16173264,3,41.513811,-82.078243,5266,258.5
16173264,7,41.515101,-82.071902,5180,258.5
16173264,11,41.516373,-82.065547,5105,258.5
16173264,15,41.517608,-82.059206,5036,258.5
16173264,19,41.518844,-82.052913,4970,256.6
16173264,23,41.520053,-82.04672,4905,252.7
16173264,27,41.521254,-82.040638,4846,246.9


texts: osm2pgrouting ERROR: relation "__ways****" does not exist

texts:
I'm using osm2pgrouting. It works for small datasets very well but when I run it on a large dataset, it start good but after some processing like 40% it starts giving following error
[*******************|                               ] (38%) Total processed: 1760000     Vertices inserted: 395 Split ways inserted 510
[*******************|                               ] (39%) Total processed: 1780000     Vertices inserted: 228 Split ways inserted 250
[*******************|                               ] (39%) Total processed: 1800000     Vertices inserted: 200 Split ways inserted 262
[********************|                              ] (40%) Total processed: 1820000
ERROR:  relation "__ways8380" does not exist

While processing FROM 1800000th          to: 1820000th way
count1820000 While processing FROM 1800000th     to: 1820000th way
[********************|                              ] (40%) Total processed: 1840000
ERROR:  relation "__ways8380" does not exist

While processing FROM 1820000th          to: 1840000th way
count1840000 While processing FROM 1820000th     to: 1840000th way
[********************|                              ] (40%) Total processed: 1860000
ERROR:  relation "__ways8380" does not exist

While processing FROM 1840000th          to: 1860000th way
count1860000 While processing FROM 1840000th     to: 1860000th way
[********************|                              ] (41%) Total processed: 1880000
ERROR:  relation "__ways8380" does not exist

While processing FROM 1860000th          to: 1880000th way
count1880000 While processing FROM 1860000th     to: 1880000th way
[********************|                              ] (41%) Total processed: 1900000

I'm running it on PostgreSQL13. Any ideas why it is giving this error.


texts: How to delete Field From Shapefile by Arcpy

texts:
I am using the following code to delete a field from a shapefile but I am getting an error!

def main():
    # Import arcpy module
   import arcpy

# Local variables:
   point_shp = "G:\\Temp\\All_Provinces.shp"

# Process: Add Field
   #arcpy.AddField_management(point_shp, "NEWFIELD", "TEXT", "", "", "25", "", "NON_NULLABLE", "NON_REQUIRED", "")

   arcpy.DeleteField_management(point_shp,NEWFIELD)
if __name__ == '__main__':
    main()

can you please let me know what I am doing wrong?


texts: How to use matplotlib with ArcGIS 10.1 without Microsoft Visual C++ Runtime Error?

texts:
I have problems using the matplotlib for python with ArcGIS 10.1

If I run the code below as a script tool in ArcMap. No pdf is generated and when I close ArcMap this error is shown.

If I run the code from IDLE a similar error is shown for C:\Python\ArcGIS10.1\pythonw.exe but the pdf is generated correct.

If I run the code in a cmd prompt no error code is shown and the pdf is generated correct.

Does anybody know how to make the script work as a script tool in ArcMap?

from matplotlib import pyplot

fig = pyplot.figure(figsize=(8.27,11.69))

axis1 = fig.add_subplot(311)
axis1.plot(range(10))

fig.savefig('multipleplots.pdf')


texts: Saving content of map.div as PNG using widget with ArcGIS Web AppBuilder?

texts:
I basically need to take a screen capture of the div WITHOUT all of widgets and panels on top of it. My plan is to export map scenes for eventual use in PowerPoint, so I want a 'what you see is what you get' look. The PrintTask() seem slow and the results always look like blurry garbage. Maybe I'm not calling it properly, but I haven't been able to get the extent the same as the screen or have the exact same resolution.

I've also tried building a canvas and adding all of the actual raw map tiles to it, but browser security won't allow me to export the 'dirty' canvas as a png due to Cross-Origin Resource Sharing (CORS).

Ideally I'll have a button in my widget called 'export' that'll generate the png and prompt the user to save it.

Anyone have any ideas?


texts: How can I open/view this set of ArcGIS files?

texts:
I have five files:

NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.img
NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.img.aux.xml
NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.img.vat.dbf
NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.img.xml
NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.lyr
NCC_Aerial_SanFrancisco_1999_2010_KelpPersistence.rrd

I understand that these files make up a raster layer. I want to be able to view them in ArcGIS. However when I try using the web developer portal, it seems to want a shapefile. This is not a shapefile. Perhaps I have to be using ArcGIS Desktop?

I can't find any documentation on this specific set of files online.


texts: Manually creating or copying an auxiliary table in QGIS 3 with Python

texts:
I have a local layer which gets updated every now and then by being copied from a WFS. The local layer has some labels shown and some hidden through auxiliary table, because QGIS seems to only accept a very specific auxiliary table to show/hide labels through labels toolbar.

What I am trying to do is to create or copy the existing auxiliary table from the old local layer to the newly created local layer with python. I've been trying some methods from here. But with little success.

I know that I can also create a Boolean field in a shapefile and point there for show/hide label. But doing that makes me unable to use show/hide label button in the toolbar. So that is not ideal.


texts: QGIS Legend - distribution of items across columns

texts:
I'm a new QGIS user, still finding my way around. I've created a map and am now preparing for export in print composer.

I've added a legend with ten items, and I want the items to break across five columns evenly. I've specified the number of columns and checked Split Layers, but the end result is three columns of two items, one column of three items, and a final column of one item:

How can I correct this to display evenly?


texts:
I need to convert raster to vector, I have data from a field in two layers. In the first i have DSM (digital surface model) and DTM (digital terain model) in the second layer. Based on difference in those
layers I see areas damaged by animals.
The program I use is SAGA.
The differences are light blue areas in attached images.
How can I convert these areas to vector?
https://i.stack.imgur.com/d46e8.jpg


texts: No result from pgRoutingLayer tsp(euclid) algorithm

texts:
I use QGIS plugin pgRoutingLayer for displaying the shortest route. Dijkstra algorithm works fine.

I tried tsp (euclid) algorithm but it drew nothing (even no error message). I thought it may be something wrong with my data. So I uploaded pgRouting sample data into my database (http://docs.pgrouting.org/dev/doc/src/developer/sampledata.html) and tried tsp on them. But result was the same (no route and no error message).

So I debugged the pgRoutingLayer plugin, caught the sql query. It looks like this and it results in 0 records

SELECT seq, id1 AS internal, id2 AS node, cost FROM pgr_tsp('
               SELECT DISTINCT id, x, y FROM
                   (SELECT DISTINCT source AS id, x1::float8 AS x, y1::float8 AS y FROM edge_table
                   UNION
                   SELECT DISTINCT target AS id, x2::float8 AS x, y2::float8 AS y FROM edge_table)
                   AS node WHERE node.id IN (7,8,11)',
               2, 12)

I may be missing some basics, I am a routing newbie. Can anyone give me a hand please? Thanks.


texts: Using String Functions (LEFT, RIGHT, MID) in Select By Attributes against a File Geodatabase using ArcMap 10

texts:
I finally (albeit reluctantly) made the switch from a personal geodatabase (pGDB) to a file geodatabse (fGDB) not too long ago.  The main reason it took so long to switch has to do with how Select by Attributes deals with string functions.

For example, I used to be able to easily select all records where the first character of a particular field starts with "A".

left([Some_Field],1 = "A")

This works without a hitch when querying against a pGDB.  However, try to do the same thing against a fGDB.

left("Some_Field",1 = "A")

This returns the following error:

There was a problem selecting
An invalid SQL statement was used.

I realize I can use the field calculator to calculate a new field equal to the first character of some other field.  But who wants to do that every time?

Is there a way to use string functions inside Select By Attributes?


texts: Eliminate small area pixel blobs of a two-value image in Google Earth Engine

texts:
I have tried to extract water body by calculated SDWI using Sentinel-1 images in Google Earth Engine :
var point = ee.Geometry.Point([-81.19416175071366,42.14552897443303]);
var sentinel_col = ee.ImageCollection("COPERNICUS/S1_GRD");

var image = sentinel_col.filterDate('2020-05-01','2020-06-01')
                        .filterBounds(point);
var sdwi = image.map(function(img){
  var wi1 =img.expression('(10*VV*VH)',
                          {'VV':img.select('VV'),'VH':img.select('VH')})
                          .rename('Water1');
  var wi2 = wi1.log();
  var wi3 = wi2.expression('water-8',{'water':wi2.select('Water1')}).gt(0.35)
                          .rename('sdwi');
  return img.addBands(wi3) .select('sdwi')
});

Map.addLayer(sdwi,[0,1],'sdwi');
Map.centerObject(sdwi);

The result seems not bad, but there are lots of small area pixel blobs（similar to noises） in the result image.
I want to eliminate these numerous pixels without losing the detail of edge (so filter kernel ways were passed).
Although there is a function called "Segmentation Image" in ENVI have achieved my thought, I am still thinking if GEE could do the same process.


texts: Correct translation of Norwegian landcover classes derived from WMS?

texts:
I am making a map of the study site located in Norway. I have accessed the awesome and free WMS data regarding landcover (I guess: N250 Arealdekkeflate) https://openwms.statkart.no/skwms1/wms.topo4.graatone?

I have made a map and added a WMS legend, but, the legend is in norwagian, of course. Please, where can i find the correct terms in English? I think my confusion originates from the the difference between land cover and land use, as well as different levels of aggregation of landcover classes by spatial scale representation. The terminology is very field-dependent and therefore simple translation by any dictionary just gives a general idea, but not the exact scientific term.

Please, what do you think about my translations? Can you suggest me a better guide?


texts:
I have a QGIS project with a large number of layers. Is it (without using Postgis) possible to define  "Views" which include different subsets of the layers?
For example:
View 1: Layers WMS topographic, roads and villages
View 2: Blank background with only village outlines and existing sewer systems
etc.
This would come in handy because I have a number of print composers which are based on different selections of the layers. If I want to make a change to one of the print composers I have to remember which layers were visible when I locked the layers in the print composer.


texts: Windowed read of a raster producing wrong projection?

texts:
I'm using rasterio to read in a .tif file and converting it to vectors using the rasterio.features module. For a particularly large .tif file, I thought I'd use the windowed read option, but for some reason the projections are coming way off when doing the windowed read. Notice how the Y coordinates are very different between the two outputs but the X coordinates are similar. Any ideas on why this is happening? See the example below:

With windowed read

with rasterio.open('raster2vector/clipped_fire2_reprojected.tif') as dataset:
    for block_index, window in dataset.block_windows(1):
        mask = dataset.read_masks(1, window=window)
        image = dataset.read(1, window=window)
        results = (
            {'properties': {'pixel_value': v}, 'geometry': s}
            for i, (s, v)
            in enumerate(
                shapes(image=image, mask=mask, transform=dataset.transform)))
        geoms = list(results)
        if not geoms:
            continue
        print(dataset.crs)
        print(geoms[0])
        break

produces

CRS({'init': 'epsg:2163'})
{'geometry': {'coordinates': [[(-1861517.9955081018, 916182.1706994411),
    (-1861517.9955081018, 915914.427799729),
    (-1860447.0239092538, 915914.427799729),
    (-1860447.0239092538, 916182.1706994411),
    (-1861517.9955081018, 916182.1706994411)]],
  'type': 'Polygon'},
 'properties': {'pixel_value': 6.0}}

Without windowed read

with rasterio.open('raster2vector/clipped_fire2_reprojected.tif') as dataset:
    mask = dataset.read_masks(1)
    image = dataset.read(1)
    results = (
        {'properties': {'pixel_value': v}, 'geometry': s}
        for i, (s, v)
        in enumerate(
            shapes(image=image, mask=mask, transform=dataset.transform)))
    print(dataset.crs)
    geoms = list(results)
    print(geoms[0])

produces

CRS({'init': 'epsg:2163'})
{'geometry': {'coordinates': [[(-1861517.9955081018, -457874.39062262117),
    (-1861517.9955081018, -458409.87642204505),
    (-1861250.2526083898, -458409.87642204505),
    (-1861250.2526083898, -458142.1335223331),
    (-1860982.509708678, -458142.1335223331),
    (-1860982.509708678, -458409.87642204505),
    (-1860714.7668089657, -458409.87642204505),
    (-1860447.0239092538, -458409.87642204505),
    (-1860447.0239092538, -457874.39062262117),
    (-1861517.9955081018, -457874.39062262117)]],
  'type': 'Polygon'},
 'properties': {'pixel_value': 6.0}}


texts: Converting values in datetime field to strings in text field using Python parser of ArcGIS for Desktop?

texts:
In the attribute table of my file geodatabase feature class I have a populated field called DateField of type date, and a text field called TextField.

I have not yet hit on a way to convert the dates in DateField to strings of dd/mm/yyyy format in the TextField.

I think this answer to Converting ArcGIS date value to string reading in mm/dd/yyyy format using ArcPy? should have given me enough clues, but when I try the following in the Field Calculator it results in the error that follows.

Parser: Python

Pre-Logic Script Code:

from datetime import datetime

def date_only(datetimeVal):
    formattedTime = datetime.strftime(datetimeVal, "%d/%m/%Y")
    return formattedTime

TextField =

date_only( !DateField! )

The error is:

Messages
Executing: CalculateField CountriesFC TextField date_only( !DateField! ) PYTHON_9.3 "from datetime import datetime\n\ndef date_only(datetimeVal):\n    formattedTime = datetime.strftime(datetimeVal, "%m/%d/%Y")\n    return formattedTime"
Start Time: Mon Feb 02 08:45:07 2015
ERROR 000539: Error running expression: date_only( u"13/02/2015 2:47:58 PM" )
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 4, in date_only
TypeError: descriptor 'strftime' requires a 'datetime.date' object but received a 'unicode'

Failed to execute (CalculateField).
Failed at Mon Feb 02 08:45:07 2015 (Elapsed Time: 0.03 seconds)

From the error I think that what I may be missing is the knowledge of how to convert the Unicode object that is being received by my date_only function into a datetime.date object.

Coincidentally, I noticed that there is currently a similar unanswered question on Stack Overflow.


texts: Do large vertex counts in PostGIS cause performance issues in WMS?

texts:
I am wondering if large geometry (geometry with large vertex count) in PostGIS might cause performance issues for GeoServer (or other product) in serving the WMS (not WFS)?

This question also ask about effect of large vertex count but more in PostGIS process.

Many thanks in advance.


texts:
I have multiple raster files that overlap each other. I am using the r.cross GRASS tool to find all of the unique combinations of values from my inputs. The tool works properly, and I get a new raster file with arbitrary values representing all of the unique combinations. However, in the Help tab of the tool under Support Files, it states that a "category file created for the output raster map layer describes the combinations of input map layer category values which generated each category." Where do I find this file? It is very important to me to know which raster layers and what values within those layers contributed to the output from the r.cross tool. The output is supposed to look something like this:

      category   category
      value      label
      ______________________________
         1       layer1(0) layer2(1)
         2       layer1(0) layer2(2)
         3       layer1(1) layer2(1)
         4       layer1(1) layer2(2)
         5       layer1(2) layer2(4)


texts: Determining which road segments get most use using ArcGIS Network Analyst?

texts:
To start, I should mention that I'm hoping to solve this using Network Analyst and (if I have to) Microsoft Access.  I am also a poor programmer - my knowledge of python is limited at best.

I have a dataset that contains somewhere around 20 000 individual origins and 130 individual destinations.  Each of the origins has a number of products (measured by volume) delivered to any number of destinations (let's say an average of 5 destinations for each origin).  Each product is shipped on its own truck to a single destination.

The problem I'm trying to solve is which roads have the most volume travel over them.

The first solution I thought of was route solver; I would have to iterate using some sort of crazy field identifier, but I'm having a hard time figuring out how to do it. I need to add an origin as a stop, then add the corresponding destination as a stop, then run... etc.  So my first problem is I'm not sure how to set that up to iterate properly (I can't picture the table to lay it out properly unless I make ~ 100 000 different tables!)

What I'm looking for is some way of determining the number of times each route is traveled - if I could link volumes to the lines somehow and then sum the volume by line segment, I could map the most used routes.

In a perfect world, each line segment along the route would link to a volume and I could use access to calculate the sums. Then import it back to ArcMap and categorize by values to display a map where the most used roads are thick and the rarely used roads are thin.

Chris W. (https://gis.stackexchange.com/users/28067/chris-w) did answer a very similar question here: Summing multiple routes using ArcGIS Network Analyst?, but I'll be honest - even though the process makes sense, the details are still escaping me.

Update:  I've had a coworker suggest that I might be looking at this the wrong way.  He suggested I look at it as a flow problem rather than a network problem.  I'll update again later if a solution presents itself.


texts: Polygonize a raster using GDAL JAVA bindings - passing raster pixel value to the resulting feature

texts:
I have a raster GeoTIFF file with values from 0 to let's say 5. When I run gdal.polygonize (I'm working with JAVA and GDAL bindings to JAVA) on it, I get features with corresponding polygons.

The problem is that the value of the raster pixels are not passed anywhere and I lose the information about which polygon belongs to which value. Is there any way to store this information anywhere?

PS. Sorry if the question is not clear. This is my first experience with GIS world (and with advanced programming also to be honest :) )


texts: How to export only one band from an image using GDAL?

texts:
I have a multi-band GeoTiff, and I want to extract only the first band and write a new image consisting of only that band.

How do I do this using GDAL?


texts: Performing transformation from GHANA GRID to WGS 84 using QGIS?

texts:
I have a shapefile which was referenced using Ghana's grid system.

I imported it to QGIS and specified WGS 84 but the coordinate bar still displays the Ghana grid values in the coordinate bar.

Is there a way to perform coordinate transformation in QGIS?


texts: Make Feature Layer from Query Layer yields an ArcPy error?

texts:
My project starts from a MXD containing Query Layers. I have to prepare the data in this "MasterMXD" for use in an offline application. Working with ArcMap 10.2.1 for this project. My Python script follows a simple structure:

Loop over all layers in the MasterMXD.
MakeFeatureLayer from the Query layer.
FeatureClassToFeatureClass to store feature layer as a feature class.
Continue from there to CreateRuntimeContent later in the script.

Here is part of my Python script:

for df in arcpy.mapping.ListDataFrames(master_mxd):
   for lyr in arcpy.mapping.ListLayers(master_mxd, data_frame = df):
      lyrName = lyr.name.encode("ascii", "ignore")
      lyrName = lyrName.replace(" ", "")
      featureLyr = arcpy.MakeFeatureLayer_management(lyr, lyrName)
      arcpy.FeatureClassToFeatureClass_conversion(featureLyr, outFolder, lyrName)

The problem I am facing is that some of the layers in the MXD pass through the script without problem, while others give an error 000229:

  ERROR 000229: cannot open GPL0 for the MakeFeatureLayer_management tool.

It is always the same layers that give the error.

I cannot see any differences between layers that give this error and the layers that do not give it. I did not construct the database that serves the query layers and cannot find any 'GPL0' anywhere. I'm puzzled!

Does anyone have an idea what is causing this ERROR 000229 showing up for some of my query layers in the MakeFeatureLayer_management tool?


texts: Missing WMS legend graphic when exporting pdf using PyQGIS

texts:
I'm developing a PyQGIS processing script to automatically load some WMS layers and export them as part of a pdf map. I import the layout template from a file, which has the QgsLayoutItemLegend already created and then add the layers I want to show in the legend. Whilst the layers load correctly and are properly rendered on the QgsLayoutItemMap in the exported pdf, the WMS legend graphic is missing.
The legend does include the WMS layer name (circled blue in the below image), just not the graphic itself. The pdf map also includes a vector layer which is correctly shown in the legend (circled red).

I've replicated this behaviour when manually entering the layers to the layout's legend item in Print Composer. However, when I adjust the size of the legend, the WMS legend graphic then appears. I tried doing this with PyQGIS in the script as a bit of a hack, but it didn't work. I've also tried  legend.updateLegend() and legend.refresh() but no luck.
Below is a code snippet of what I'm doing...
#read layout template
layoutFile = open('layout_path')
layoutContent = layoutFile.read()
layoutFile.close()
document = QDomDocument()
document.setContent(layoutContent)

currentProject = context.project()
layout = core.QgsPrintLayout(currentProject)
layout.initializeDefaults()
layout.loadFromTemplate(document, core.QgsReadWriteContext(), False)

#get legend object and add layers to it
legend = layout.itemById(legend_id)
legendModel = legend.model()
legendRootGroup = legendModel.rootGroup()
legendRootGroup.clear()

for l in layers_to_show: #a list of QgsVectorLayers and QgsRasterLayers
    legendRootGroup.addLayer(l)

legend.refresh()
legend.updateLegend()

#export pdf
exporter = core.QgsLayoutExporter(layout)
exporter.exportToPdf(saveAs, core.QgsLayoutExporter.PdfExportSettings())


texts: R and Mapview: How can I chose different colors for different layers?

texts:
Here is an example:

library(sp)
library(mapview)

data(meuse)

d<-st_as_sf(meuse,coords =  c("x","y"))
d1<-d[sample(nrow(d),20),]
d2<-d[sample(nrow(d),20),]
mapview(list(d1,d2))

Both layers, d1 and d2, and displayed using the same color. How can I chose different colors, say red for d1 and blue for d2?


texts: Google Earth Engine Sum Area of Image and time series

texts:
I'm calculating NDSI from landsat-8 image and I want to show the snowy area through a chart. I'm able to calculate the area for each pixel based on NDSI but I don't know how to sum the pixel values... here's the code:
//CALCULATING NDSI FUNCTION
function addNDSI(image){
  var NIR = image.select('B5').gte(0.1).rename('NIR');
  var RED = image.select('B4').gt(0.12).rename('RED');

  var NDSI = image.normalizedDifference(['B3', 'B6']).rename('NDSI');
  var NDSImask = NDSI.gt(0.3).and(NIR).and(RED).rename('NDSImask');

  NDSI = NDSI.multiply(NDSImask);

  //CALCULATING SNOWY AREA in KM^2
  var areaImage = ee.Image.pixelArea().divide(1000000); //from meters to kilometers (^2)
  var snowAreaImage = NDSImask.multiply(areaImage).rename('snowArea');

  return image.addBands(NDSI).addBands(snowAreaImage);
}

//SELECTING AREA OF INTEREST & FILTERING: PATAGONIA
var regionPatagoniaFiltered = L8.filterBounds(roiPatagonia)
    .filterMetadata('CLOUD_COVER','less_than',40)
    .filterDate('2014-01-01','2018-12-31');

//DISPLAYING SATELLITE IMAGES OF THE INTEREST REGION ON MAP
var medianPatagoniaImages = regionPatagoniaFiltered.median();
Map.addLayer(medianPatagoniaImages,mapImageParams,"Layer - Average RGB image in this period");

var with_NDSI = regionPatagoniaFiltered.map(addNDSI);

Map.addLayer(with_NDSI.median(),mapNDSIparams,"Layer - Average NDSI in this period ");
print(Chart.image.series(with_NDSI.select('NDSI'),roiPatagonia));

var chart = ui.Chart.image.series({
  imageCollection: with_NDSI.select('snowArea'),
  reducer: ee.Reducer.sum(),
  region: roiPatagonia,
  scale: 30,
});
print(chart);

How you can see, I tried to use the reducer.sum, but it seems to not work...


texts: Using Label Expression with multiple fields and assign different colors

texts:
I'm trying to use an expression to label a layer by using multiple fields, and assign different colors using QGIS 2.16.
In ArcGIS I am able to achieve this using the following expression:
def FindLabel([DesceasedName], [ownerName], [all_spaceid]):

  if [DesceasedName] = "":
    return "<BOL><FNT size = '5'>" + [DesceasedName] + "</FNT></BOL>"+   '\n' + "<CLR red='115' green='38' blue='0'><BOL>" + [all_spaceid] + "</BOL></CLR>"
  elif [DesceasedName] == "":
    return "<CLR blue='255'><FNT size = '5'>" + [ownerName] + "</FNT></CLR>"  + '\n' + "<CLR red='115' green='38' blue='0'><BOL>" + [all_spaceid] + "</BOL></CLR>"
  else:
    return ""

How to achieve the same labeling in QGIS?


texts: QGIS, Postgis: Is it possible to create a new postgis layer?

texts:
Situation:

I have a postgis vector lay that has lots of line segments (the features).  I want to create a set of new layers each of which contains a subset of the features of the original layer and have these layers also stored back in postgis.  Ideally with the actual features stored just once.

The only way I can see to create a new layer as opposed to load one from some source is via Layer->new which gives an option of shapefile or spatialite.

Choosing shape file I first tried to specify gid and fid attributes for the new layer but SPIT complains about duplicate field gid when I try to load it into postgis.

If I do not create any attributes for the new layer SPIT happily loads it, presumably duplicating all the features.

I am guessing this is a 'feature' rather than an oversight because if it were possible to share features between layers then then what would happen if a feature were split? If things were to work properly then all layers with referenced that feature would need updating with the new gid of the new portion.

Is this interpretation correct?

Russell


texts:
I need to examine certain areas over a period of several years to view the effects of urbanisation. It's largely a visual exercise. I've found one tool that does this and is easy to use[1] but it only provides the historic images for the most zoomed-out levels.

Can anyone recommend a resource to achieve this? We have no budget, but I'm open to paid services if no free ones are available.

PS: If there is specific terminology I should be aware of to help with internet searches, please let me know. The important thing is the record of historical images to compare land changes over time such as deforestation. The location of interest this time is Australia.

1: https://zoom.earth/#-33.340544,146.887772,8z,sat,am,2015-05-06


texts: GEOS segmentation fault when loading libc.so.6

texts:
I'm trying to run a program that uses shapely on a GCE Virtual Machine (Ubuntu 18.04). I'm getting an error on the geos.py file, here is the relevant trace I'm getting:

[geos:73] - [DEBUG]: Found GEOS DLL: <CDLL '/home/loic/venv/lib/python3.6/site-packages/shapely/.libs/libgeos_c-5031f9ac.so.1.13.1', handle 31283c0 at 0x7f5aefc58940>, using it.
[geos:38] - [DEBUG]: Trying `CDLL(libc.so.6)`
Segmentation fault (core dumped)

I'm running this program inside a python virtualenv, with Shapely==1.7.0. The thing is, when I'm launching this program locally, with basically the same environment, I'm not getting any error, CDLL(libc.so.6) is working fine.

Here's the output of ls -l /lib/x86_64-linux-gnu/ | grep libc:

-rwxr-xr-x 1 root root 2030544 Apr 16  2018 libc-2.27.so
lrwxrwxrwx 1 root root      12 Apr 16  2018 libc.so.6 -> libc-2.27.so
lrwxrwxrwx 1 root root      18 Sep 24  2017 libcap-ng.so.0 -> libcap-ng.so.0.0.0
-rw-r--r-- 1 root root   18712 Sep 24  2017 libcap-ng.so.0.0.0
lrwxrwxrwx 1 root root      14 Nov 16  2017 libcap.so.2 -> libcap.so.2.25
-rw-r--r-- 1 root root   22768 Nov 16  2017 libcap.so.2.25
-rw-r--r-- 1 root root  190784 Apr 16  2018 libcidn-2.27.so
lrwxrwxrwx 1 root root      15 Apr 16  2018 libcidn.so.1 -> libcidn-2.27.so
lrwxrwxrwx 1 root root      17 Jan 22 14:40 libcom_err.so.2 -> libcom_err.so.2.1
-rw-r--r-- 1 root root   14248 Jan 22 14:40 libcom_err.so.2.1
-rw-r--r-- 1 root root   39208 Apr 16  2018 libcrypt-2.27.so
lrwxrwxrwx 1 root root      16 Apr 16  2018 libcrypt.so.1 -> libcrypt-2.27.so
lrwxrwxrwx 1 root root      23 Aug 23  2018 libcryptsetup.so.12 -> libcryptsetup.so.12.2.0
-rw-r--r-- 1 root root  310040 Aug 23  2018 libcryptsetup.so.12.2.0

Has anyone experienced this before, or know how to fix this issue? I can provide more information if needed.


texts: Extracting data using arcpy.da.SearchCursor?

texts:
I have a shapefile of point features showing fire locations in a given year. One of the fields is named "YYYYMMDD" and indicates the corresponding date that that fire happened. In Python, I am trying to use the arcpy.da.SearchCursor function to append each date into a list and find the most recent fire by date.

My code looks like this:

listname = []
cursor = arcpy.da.SearchCursor(fires, "YYYYMMDD")

for row in cursor:
    listname.append(row)

print 'the first fire date is ',min(listname),'.'
print 'the last fire date is ',max(listname),'.'

What prints for min and max listname is: (20010515.0,) and (20011002.0,)
But what I would like to get as an output is: 20010514 and 20011002

How do can I get that as my output?

I tried indexing what I got to exclude what I don't need but that doesn't work.


texts: Automate CSV to Point Layer to Polygon (Voronoi) on QGIS (Mac or Windows)

texts:
I have a map I update on a monthly basis, following these steps:

Import new CSV with lat/long and other values as point layer
Generate Voronoi polygons for the point layer
Clip the Voronoi polygons to the base layer
Dissolve the clipped polygon layer based on values in Column C
Save the resulting polygon layer for use as a base map in other projects

Since the steps are the same every time, I am sure there is some way to automate this, but I am weak on writing a script to run QGIS from the command line. I primarily use Mac, but I also have access to a PC for this project. Currently running QGIS 2.18.2

My question is similar to the unanswered Automate the CSV importing process so that map is updated dynamically in QGIS?, so I hope I have provided more specific details than the author of that question to get a clear answer.


texts: gdal_merge shifts rasters after processing, layer stack is slightly off

texts:
I have experienced that gdal shifts rasters. The outputs of gdal_warp are often skewed, have different pixel sizes and so forth. A known error which has also been discussed here.

But now, I had to realize, that gdal_merge also shifts rasters! Well it is only some mm of whole MODIS tiles (pixel size 250m), but I can't imagine what is the reason for it, all layers have the same extent, the same coordinate system and so forth. Why should gdal_merge shift inputlayers and not simply merge or stack them?

[code example:

gdal_merge -separate -of GTiff -o D:/MODIS_TIF_NDVI_RECTCLIP_YEARLY_STACK/2005_PARAMETERS_NDVI_STACK.tif D:/MODIS_TIF_NDVI_RECTCLIP_SEASONAL_PARAMETERS/2005_SOS_s1 .....]


texts: Incremental retrieval of data from a database - OpenLayers

texts:
I have a new task now:

The user will input the lat-lon values in a HTML form and when an INSERT button is clicked these values should be stored the database. Now when a submit button is clicked the markers of all the lat-lon values stored in a table in the database should be visible on the map. This is fairly not a tough task, but... when the user provides next value of a point, this should be stored in the same table and this particular point marker should be added to the map (the map already contains markers of previous points)... that is, instead of reading the whole data (including the new point) from the database and projecting on the map, since before adding new point to the database we have already pointed the markers of previous values, now only the new point marker should be added to the map.

For better understanding here is the example...

I already have 999 points in my database and these are pointed on the map. Now the user gives 1000th point values and this will be stored in the DB. now if we reload (click on submit button) the page, instead of reading all the 1000 points from the database and create markers, it should only read the new point and add this to already existing map with 999 markers.

I have a simple logic for this, table1 in the database contains all the values and table2 only the new value. Whenever the user insert a points it is simultaneously added in both the tables and when the user submits the only value from table2 is deleted (so that next new value can be stored). So reading just from table2 works here, but how to add this point to the already existing map is the main task.

I don't have any idea whether this can be done in OpenLayers or not, so please help me how to do this.

Hope my intention is understandable.

Till now I have not worked with databases in OpenLayers, so  please suggest me which database (like MySQL, etc...) suits better for this task... and any other suggestions on storing/retrieving to/from the database.


texts: Identifying geographies that match certain climatic profiles through GIS

texts:
The long and short of this research is trying to find roads that lie within geographic areas that fit a given climatic profile. In this case, I need to find areas within the United States with relatively cold falls, winters, and springs that get snow. These places may not have a lot of snowfall on an annual basis, but the snow persists for longer periods due to the colder weather patterns.

Are there any websites that already collect similar data in a format ready for GIS analysis? If not, what should I look for when gathering data? There may not be a singly perfect answer for this question, but I mostly am seeking direction to refine my search appropriately.


texts: Choosing between several very similar projections

texts:
When working with data from different sources, they often come in with different EPSG codes, even when they are using the same nominal projection.  For example, I often work in Washington State Plane North NAD83 in US feet.  Using QGIS I generally select EPSG 2285, but when I'm working with a shapefile generated in ArcMap, it often comes in on EPSG 102748.  I recently opened a geodatabase also in Washington State Plane North NAD83 in US Feet that was in EPSG 6597.  A search on EPSG.io also lists EPSG 32148 as a Washington State Plane North NAD83 in US Feet.

When using QGIS all of these appear to perform the same, and I haven't run into issues with with strange results because of different projections.  None of these projections appear to be deprecated, and I haven't found anything to recommend one vs another.  Is there any reason to prefer or avoid any of these, or am I getting concerned over trivia?


texts: How to interpolate subsurface geology from borehole data

texts:
I have borehole data (x, y, z values for each subsurface measurement, and its associated lithology). I would like to interpolate the borehole data to create 1) interpolated subsurfaces and 2) lithological profiles at specific cross-sections. I only have access to the ArcGiS package w/ extensions, but no other software packages, nor the budget to access them.

I need to know if the above is possible with the tools and data that I have, and if so, what tools and methodology I should be looking into.


texts: Expanding extent of file geodatabase feature class in ArcGIS Pro?

texts:
An answer by @CurtisPrice to Recalculating full extent of File Geodatabase feature class? alerted me to the Recalculate Feature Class Extent in ArcGIS Desktop which is also available in ArcGIS Pro.

Using ArcGIS Pro 1.4.1, and a file geodatabase feature class, I would like to be able to do something similar to this tool, which reduces the feature class extent, but in my case I would like to increase the feature class extent.

The reason I want to do this is because while using ArcPy to move point features via arcpy.da.UpdateCursor and their SHAPE@XY tokens I am triggering this error:

  RuntimeError: The coordinates or measures are out of bounds.

This is because the location that I am moving the point to is not only west of its original location but also west of any other points in the feature class.

Is there a tool which can be used to expand the extent of a feature class?

If there is not, a workaround that I will try will be to create a feature class of four points that are all outside of the current extent, Append them into the feature class to hopefully trigger an update of its extent, do my point movements, and then Select all but those four points out into a new moved feature class.


texts: Which projection is best for mapping the contiguous United States?

texts: If I wanted to project latitude & longitude data for the contiguous United States (the United States excluding Alaska and Hawaii), which projection would I use? I prefer more accurate distances followed by shapes


texts: Check if points that share attribute intersect polygons

texts:
How to check whether points started in red and ended in a blue polygon (and vice versa) while sharing the same attribute?

To check if the first location of the unique ID is in one polygon and last in another.

Points sample:


texts: Show/hide markers depending on zoom-level

texts:
I would like to show icons if the user zooms to level 7 or more and hide the markers when zoom level is 6 or lower. I've tried using the 'zoomend' property in my function like so:

var shelter1 = L.marker([55.1, 11.6], {icon: shelterIcon}).addTo(map);

    map.on('zoomend', function() {
        var currentZoom = map.getZoom();
        if(currentZoom >= 7) {
          shelter1.bindPopup("Shelter");
        }
        else {
          shelter1.hide();
        }
    });

I know that hide() is not valid. But I need something similar. Couldn't find anything in the Leaflet documentation.


texts: issue with WEB ADF based Web application

texts:
I have a WEB ADF based application that I am having some issues with. When I start app in debug mode, it shows me a blank page but I can right click and choose the option view in  browser and it works fine and loads properly. I am trying to figure out what I need to do to load this page propely in a debug mode. Appreciate any ideas.

Thanks
Jay


texts: How to change values in a raster file at particular positions?

texts:
I have a collection of coordinates which comes from a QGSGeometry object. The shape is a polygon.

example:

qgsGeometry [[(75.2172,17.6566), (75.8499,17.9613), (75.5218,17.2113), (75.2172,17.4574), (75.2172,17.5277), (75.2172,17.6566)]]

How can I change the values of a raster file (geotif) only at the raster cells inside of the polygon (red line in screenshot)?

I know how to transform a raster and write or update with gdal it. But I don't know how to change values only at an given area like the inside of a polygon.


texts:
Is there a KML tag or specification by which you can force a given layer (specifically, a screen overlay) to be "always on", in that there is no way for the user viewing the file to make that layer invisible, so long as they are viewing anything in that file? I'm looking in particular for a solution that will work when Google Earth is used to view the KML file.

That is, if my KML file has multiple features and layers, one of which is called "Overlay", say, then if the user loads my KML file and turns on any of the layers in it, then Overlay will automatically turn on. If the user tries to turn Overlay off while something else in the file is still visible then either (a) Overlay cannot be turned off, or (b) everything else in the file is made invisible along with Overlay.

Of course I have no problem with the user removing my KML file from the GE session, or disabling it entirely, but if any of it is shown, then Overlay should be as well.


texts: Openlayers Vector Layers with Google Map Base Layer (via Grid): Vectors Stick but can't get correct resolution/projection?

texts:
Revamping this question instead of opening a new one.

Previous issue: vector layers were not "sticking" to the map, so when you would zoom in or out the vector layers would "come out" at you. Note their lat and long didn't change, but it looked pretty tacky when zooming-as they would not stick to the map.

It seemed like the only way to do this (the stick factor) was to use a WMS or Grid Layer: correct me if I am wrong. I chose Grid in this instance.

Note: the "stick factor" is not an issue when the base layer is OSM or OpenLayers, only Google base layers. Here is an example demonstrating this at bottom of page:

The code I currently have for this Grid is currently this:

var options = {
            resolutions: [1.40625,0.703125,0.3515625,0.17578125,0.087890625,0.0439453125],
            maxExtent: new OpenLayers.Bounds(-180, -90, 180, 90),
            //restrictedExtent: new OpenLayers.Bounds(-360, -90, 360, 90),
            numZoomLevels: 8,
            singleTile: true,
            isBaseLayer: true,
            wrapDateLine: true,
            controls: [
            new OpenLayers.Control.OverviewMap(),
            new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.ScaleLine( { maxWidth: 150 })
            ],
            getURL: function() {
                    var center = this.map.getCenter(),
                    size = this.map.getSize();
                    return [
                    this.url, "&center=", center.lat, ",", center.lon,
                    "&zoom=", this.map.getZoom(), "&size=", size.w, "x", size.h
                    ].join("");
            }
    };

Now the first thing that comes to mind is that I feel I should have an options variable for the map object and an options variable for the Grid layer (what are your thoughts on this?), as shown below:

 var grid = new OpenLayers.Layer.Grid(
        "Google Streets",
        "http://maps.googleapis.com/maps/api/staticmap?sensor=false&maptype=roadmap",
        null, options
    );
    map.addLayers([grid]);

The left picture below is how the resolution was before: just using a standard OpenLayers Google layer.

The picture to the right is using the Grid layer and point to the Google default layer: mentioned in the code above.

I just posted a snippet of Africa to show the resolution issue, and as for the points (all my points are in the US, they are completely out of whack also (the left most point on the map should be on the west coast of the USA, while the right most point should be on the east coast of the USA.:

1) Should I create a separate options variable for the grid layer?

2) Is adjusting the resolution within the options the correct way to fix this?

3) As you can see in getURL there is no mention of projects or transformations: how do I account for the points showing up correctly in the USA?


texts: Get path of feature classes in feature datasets without arcpy

texts:
I'm looking for a way to get the full path of a feature class within feature datasets (either in personal or file geodatabases) without using arcpy.
My first attempt was to use ogr which works fine if the feature classes are stored in the geodatabase directly. But I didn't manage to get the full path if a feature dataset is involved. I tried the "OpenFileGDB" driver as well as the "FileGDB" driver.
My test geodatabase stores five feature classes, four of them are within the feature dataset and one directly in the geodatabase.

My code so far:

gdb_path = r"C:\my_path\test_gdb_fd.gdb"
gdb_driver = ogr.GetDriverByName("FileGDB")
gdb = gdb_driver.Open(gdb_path)
for featsClass_idx in range(gdb.GetLayerCount()):
    featsClass = gdb.GetLayerByIndex(featsClass_idx)
    print featsClass.GetName()

this prints the basenames of all feature classes without any information whether the feature class is within the feature dataset or not.
I'm using GDAL 2.2.4, FileGDB API 1.3, Python 2.7 32 bit


texts: How to turn off geojson layers in Leaflet when the site starts?

texts:
I have a problem with my geojson layers. When I lunch my site they are all switched on, but I don't want to. Is there any solution to have lunched only one layer when my site starts, not the all one? Here is the simple example what I want https://jsfiddle.net/expedio/dtLyvffg/

And this is the code with my layers:

var controlLayers = L.control.layers().addTo(map);

$.getJSON("data/Zabudowa.geoJSON", function(geojson) {
  var geojsonLayer = L.geoJson(geojson, {
    style: function(feature) {
      return {
        'weight': 0,
        'fillColor': 'brown',
        'fillOpacity': 1
      }
    },
    onEachFeature: function(feature, layer) {
      var popupText = "Numer obiektu: " + feature.properties.OBJECTID + "
Sposob uzytkowania: " + feature.properties.type + "
More info";
      layer.bindPopup(popupText);
    }
  }).addTo(map);
  controlLayers.addOverlay(geojsonLayer, 'Zabudowa');
});

controlLayers.addOverlay(geojsonLayer, 'Drogi');
});

var myIcon = L.Icon.extend({
  iconUrl: 'icons/cemetery-24.png'

});
$.getJSON("data/miejsca_pamieci.geoJSON", function(geojson) {
  var geojsonLayer = L.geoJson(geojson, {
    pointToLayer: function(feature, latLng) {
      return new L.Marker(latLng, {
        icon: new myIcon({
          iconUrl: 'icons/cemetery-24.png',
          iconSize: [24, 24],
          iconAnchor: [12, 24],
          popupAnchor: [0, -24]
        })
      })
    },
    onEachFeature: function(feature, layer) {
      layer.bindPopup("Opis: " + feature.properties.Opis)
    }
  }).addTo(map);
  controlLayers.addOverlay(geojsonLayer, 'Miejsca_pamieci');
});

var myIcon = L.Icon.extend({
  iconUrl: 'icons/accessdenied.png'

});

$.getJSON("data/Fotografie.geoJSON", function(geojson) {
  var geojsonLayer = L.geoJson(geojson, {
    pointToLayer: function(feature, latLng) {
      return new L.Marker(latLng, {
        icon: new myIcon({
          iconUrl: 'icons/accessdenied.png',
          iconSize: [32, 37],
          iconAnchor: [16, 37],
          popupAnchor: [0, -37]
        })
      })
    },
    onEachFeature: function(feature, layer) {
      layer.bindPopup(feature.properties.Name, {
        maxWidth: "auto"
      })
    }
  }).addTo(map);
  controlLayers.addOverlay(geojsonLayer, 'Fotografie');
});


texts:
I am having trouble with CopyFeatures in my script. The files went being copied to the output folder. Can someone please assist. Here's a snippet of my script for ArcGIS 10.3;

#Loop through shapfiles in folder and reproject
for fc in fcList:
fcspatialRef = arcpy.Describe(fc).spatialReference.name
if fcspatialRef != spatialRef==True:
    arcpy.Project_management(fc, outFolder + "\\" + fc, template)
else:
    fcspatialRef != spatialRef==False:
        arcpy.CopyFeatures_management(fc, outFolder)

# Print shapefile Project results
        arcpy.AddMessage(fc)

# Return any errors
except:
arcpy.AddMessage(arcpy.GetMessages())


texts: Get centroid of SimpleFeatures in GeoTools8

texts:
I've got some SimpleFeatureCollections and now I want to get the Centroid of every SimpleFeature in the collection.
Unfortunally I've only found a centroid method for Geometry.

Is there any way to get the centroid?


texts: Creating Drive Time Polygons Across State Boundaries Using OSM?

texts:
I am using pgRouting on OpenStreetMaps to create drive time polygons. The challenge is that at state boundaries, these polygons effectively stop, since the roads in OSM are segmented at the state level.

Does anyone have a recommended solution for overcoming this challenge?


texts: Does raster intersect polygon in ModelBuilder?

texts: In ArcMap, I'm building a model which is iterating over ArcGrid rasters. I only want the steps to run where the raster intersects with a certain polygon. How do I add this condition to ModelBuilder?


texts: Why does Landsat 8 panchromatic band NOT include the infrared?

texts:
The panchromatic band on Landsat 8 OLI sensor (band 8) covers a narrower portion of the spectrum than the corresponding band in Landsat 7 ETM+: 0.50 - 0.68 micrometers (L8 OLI) in contrast to 0.52 - 0.90 micrometers (L7 ETM+) (see complete tables here). Most remarkably, L8 OLI band 8 does not cover the near infrared, thus precluding pan sharpening in that portion of the spectrum and making up a significant difference with the L7 ETM+ panchromatic band. Why?

In this interesting blog post, the author mentions the issue without being able to find the logic behind it, and citing the explanation given by Irons et al (2012) “The OLI panchromatic band, band 8, is also narrower relative to the ETM + panchromatic band to create greater contrast between vegetated areas and surfaces without vegetation in panchromatic images" as not clear or insufficient, and hinting to possible design glitches.

I am also left wondering about this and could not find a decent explanation anywhere. Does anyone in this forum know why?

Irons, J.R. et al. (2012) The next Landsat satellite: The Landsat Data Continuity Mission. Remote Sens. Envion.,122, 11-21. doi:10.1016/j.rse.2011.08.026


texts:
I'm working on a project of water supply where I have to merge line features representing pipes  if  they have the  same material of construction and if they are touching each other. The merge is done two by two what it means that some features will be duplicated in some cases like described in the figure below :

this shows exactly my issue. After the merge I will get three records but what I want is just one record which encompasses the whole pipes that fill the conditions put in the where clause :

Here is the query that helps me do the merge :

drop table if exists touches_material;
create table touches_material as
select   distinct a.*,
st_Union(a.geom,b.geom) as fusion  from pipe a, pipe b
where  a.id < b.id and a.material = b.material  and
st_touches(a.geom,b.geom)
group by a.id,a.geom,b.geom

the following picture picture shows the expected result on a test data, it's realized via QGIS GIS software :

but this is what I' getting with my query :


texts: Differences between triangulation and trilateration?

texts:
Looking around I noticed that many people interchange the terms (triangulation and trilateration) for the same sense.

What is the correct sense of Triangulation and what are the differences from Trilateration?