texts: Identify points in a shapefile that have no other point within x km with QGIS

texts:
This is probably a naive question but I am struggling as a new user to QGIS.

I have a very large shapefile (275,000 points, but can break this into about 10 subregions if necessary for faster processing).

I want to identify all points that have no other point within 200 metres and then code each of those points with the value "unique" in a field of the file.

For all the other points that are part of local clusters I then want to code those as "clustered".

Having achieved that, I want to then select just one for each cluster on a random basis to retain in the data set, discarding the others.

Currently I am failing to achieve step 1 so any assistance would be welcome.


texts: WFS Mapguide Maestro. Can't query features?

texts:
I'm using MapGuide Maestro 5.0 and Geoserver 2.4.0.

I have added a new WFS feature source to the viewer, and used the following url:

 http://localhost:8080/geoserver/scotwat/wfs?

In the vector layer def it shows the fields and these have been checked to be visible.

The viewer shows the features, but when I click on them to select, nothing happens. I have looked through the Properties of both the feature and layer, and cannot find an option to turn it on. I have also ensured that the box 'Features on the layer are selectable (if visible)' is checked. I also added the wfs to qgis, and can select features through there.

EDIT: I have progressed slightly, but am still stuck. The problem lies with the primary key of the WFS service. Within geoserver the option 'Expose primary key' is checked. MapGuide produces the following warning the feature source is validated:

 No primary key defined for class: scotwat:paris_hetero, features will not be selectable and feature joins on this class may produce unexpected results

The data source the WFS is generated has a primary key, but this has not been seen by MapGuide Maestro. Any help? If I know what the primary key how can I define it in MapGuide Maestro?


texts: Error while exporting data through Mappia Publisher in QGIS 3.14

texts:
I installed Mappia Publisher in QGIS 3.14 and try to share some data. It works until almost the end. It gives me a Python error. I do not know why. Also, it takes me almost two hours to run.
Finished map tile generation. Uploading changes to Github.
Github found commiting to your account.
Git: Pulling remote repository current state.
Elapsed 2s on this step.
Git: Add all generated tiles to your repository.
Elapsed 30s: Please wait, identifying changes on your repository.
Elapsed 61s: Please wait, identifying changes on your repository.
Elapsed 92s: Please wait, identifying changes on your repository.
Elapsed 114s on this step.
Git: Committing changes.
Elapsed 9s on this step.
Git: Pushing modifications to remote repository.
Elapsed 0s on this step.
Error, publication failed! Reporting error to Mappia development team.
Traceback (most recent call last):
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\mappia_publisher_algorithm.py", line 792, in processAlgorithm
return self.generate(writer, parameters, context, feedback)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\mappia_publisher_algorithm.py", line 619, in generate
GitHub.publishTilesToGitHub(writer.folder, ghUser, ghRepository, feedback, self.version, ghPassword)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\GitHub.py", line 385, in publishTilesToGitHub
GitHub.pushChanges(repo, user, repository, password, feedback)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\GitHub.py", line 353, in pushChanges
return GitInteractive.pushChanges(repo, user, repository, password, feedback)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\GitHub.py", line 682, in pushChanges
GitHub.getGitPassUrl(user, repository, password), "master:refs/heads/master")
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\UTILS.py", line 144, in runLongTask
return job.result()
File "C:\OSGEO4~1\apps\Python37\lib\concurrent\futures\_base.py", line 425, in result
return self.__get_result()
File "C:\OSGEO4~1\apps\Python37\lib\concurrent\futures\_base.py", line 384, in __get_result
raise self._exception
File "C:\OSGEO4~1\apps\Python37\lib\concurrent\futures\thread.py", line 57, in run
result = self.fn(*self.args, **self.kwargs)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\git\cmd.py", line 542, in <lambda>
return lambda *args, **kwargs: self._call_process(name, *args, **kwargs)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\git\cmd.py", line 1005, in _call_process
return self.execute(call, **exec_kwargs)
File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\mappia_publisher\git\cmd.py", line 822, in execute
raise GitCommandError(command, status, stderr_value, stdout_value)
git.exc.GitCommandError: Cmd('C:\Program Files\Git\bin\git.exe') failed due to: exit code(1)
cmdline: C:\Program Files\Git\bin\git.exe push https://liaduarte:0a18394b1e7b22cc57d853d1f920b6032dbc93e9@github.com/liaduarte/Mappia_Example/ master:refs/heads/master
stderr: 'error: src refspec master does not match any.
error: failed to push some refs to 'https://...@github.com/.../Mappia_Example/''

Execution failed after 1162.77 seconds


texts:
I'm trying to check the "fiberEquipment" table against the "material_cd" table and find instances where there are discrepancies.
E.G --
Checking that "material_cd" 10785823 in "fiberEquipment" has all the correct values associated with it (inputportcount, outputportcount, inputblockcount, etc.) per the "material_cd" table.

The code I'm using below takes the values from the "fiberEquipment" table, inserts them into a dictionary and than uses a search cursor invoking the dictionary elements to check for discrepancies. The code currently works without any errors, but the logic seems to be flawed, as it's not doing anything when all fields match.

# Fiber Equipment/Material Catalog fields
matcd_fields = ["material_cd", "inputportcount"]
fibequip_fields = ["material_cd", "inputportcount", "outputportcount", "inputblockcount", "outputblockcount", "objectid"]

# Compare FIBEREQUIPMENT table to Material Catalog table and determine which objectid's have records that aren't matching
 ipc = {} # 'material_cd' : 'inputportcount'
 ipcOID = {} # 'material_cd' : 'objectid'
 opc = {} # 'material_cd' : 'outputportcount'
 opcOID = {} # 'material_cd' : 'objectid'
 ibc = {} # 'material_cd' : 'inputblockcount'
 ibcOID = {} # 'material_cd' : 'objectid'
 obc = {} # 'material_cd' : 'outputblockcount'
 obcOID = {} # 'material_cd' : 'objectid'

with arcpy.da.SearchCursor(fiberEquipment, fibequip_fields) as cursor:
    for row in cursor:
        ipc[str(row[0])] = row[1]
        ipcOID[str(row[0])] = row[5]
        opc[str(row[0])] = row[2]
        ipcOID[str(row[0])] = row[5]
        ibc[str(row[0])] = row[3]
        ibcOID[str(row[0])] = row[5]
        obc[str(row[0])] = row[4]
        obcOID[str(row[0])] = row[5]

objectID = []

with arcpy.da.SearchCursor(material_catalog, matcd_fields) as cursor:
    for row in cursor:
        if ipc[str(row[1])] == row[1]:
             objectID.append(ipcOID[row[0]])


texts:
I'm trying to make a precision topology map in QGIS 3.14 where my map's corners are defined by geographic coords: this means it should have the same range of degrees on top as on the bottom, and same on left as on the right side of map. For example:

I use the UTM coordinate system, U34 zone, EPSG: 32634.
Coords of my corners:
Top-left: '549108, 5506822'
Top-right: '553523, 5506822'
Bottom-left: '549108, 5502204'
Bottom-right: '553523, 5502204'
Actually I feel like my map in the print-composer is oriented in this way (red square is my map range):

And my desired output should look like this (along the parallels):


texts: OpenLayers 3 Geoserver GeoWebCache WMS tile alignment issue

texts:
I've started moving to OpenLayers 3 - very nice. One issue I have is with the dynamic tile size sent to my Geoserver - which is often 333.

This doesn't line up with the cache (which is 256 for EPSG:900913) so the Geoserver has to build a new tile - somewhat slow ...

Sometimes I get 256 (e.g. from Firefox on Linux) but usually it's off (e.g. Chrome on Linux, Android, Safari on iOS).

I've looked at a number of questions on this site, and have gleaned that I might need to set a tileGrid, but how do I do this?

Do I need to set resolutions on my map?

Here's my code:

var wms =   new ol.layer.Tile({
  source: new ol.source.TileWMS({
    url: '/geoserver/wms',
    params: {'LAYERS': 'test:wms_layer', 'TILED': true, },
    serverType: 'geoserver'
  })
});

var osm = new ol.layer.Tile({
  source: new ol.source.OSM()
});

var map = new ol.Map({
    target: 'map',
    layers: [
      osm,
      wms,
    ],
    view: new ol.View({
      center: ol.proj.transform([11.953, 57.757,], 'EPSG:4326', 'EPSG:3857'),
      zoom: 16,
    })
});


texts: Adding two values into a categorized symbology in order to toggle both values on and off

texts:
I'm representing an archaeological site in QGIS and I have expressed the multiple contexts (stratigraphic units) as polygons with categorized symbology so that the client can toggle off and on the contexts they want to see etc.

My issue is that, oftentimes, these contexts overlap or cut into other contexts due to historic construction phases and deposition/accumulation.
So, I've added an attribute called 'phases' and I want to include that alongside my 'context' symbology so that if you toggle off one phase, all the contexts within that phase disappear while the phase that is visible, the individual contexts can be toggled on and off.

Is there a way for this to be done?


texts: How to create new layer within a DXF file?

texts:
Using OGR in Python, I'd like to create some geometries within a DXF file, but not in the "0" layer, but in new ones.

I tried with:

import ogr, os, osr, sys
ds = ogr.GetDriverByName('DXF').CreateDataSource('output.dxf' )
myLayer = ds.CreateLayer('entities')
newFeature = ogr.Feature(feature_def = myLayer.GetLayerDefn())
newFeature.SetGeometryDirectly(ogr.CreateGeometryFromWkt( 'LINESTRING(0 0, 100 100)'))
myLayer.CreateFeature( newFeature )
newFeature.Destroy()
ds.Destroy()
ds = None

But I don't have any idea about how to continue and create new layers.

I tried also with

layers = ds.ExecuteSQL("INSERT INTO entities (Layer) VALUES ('Layer1')")

but it does not work.

Any suggestion?
Is it worth using another DXF+Python library?

Thank you very much.


texts: Provide credentials for Oracle database when opening MXD file?

texts:
I need  to open a MXD file from a Python toolbox (so that I can automatically export it to a PDF). The file contains a connection to an Oracle 11g database, and I need to somehow provide the credentials to the database. None of the methods I have tried so far have worked - the PDF always comes out without any of the data from the database. How can I tell ArcGIS what username and password to use?

Attempt 1: Reference a SDE file:

arcpy.env.workspace = "credentials.sde"
mxd = arcpy.mapping.MapDocument(mxd_file)
arcpy.mapping.ExportToPDF(mxd, pdf_file)

Attempt 2: Create a SDE file:

arcpy.CreateArcSDEConnectionFile_management (
   out_folder_path = "C:\\Path",
   out_name = "temp.sde",
   server = "prod.world",
   service = "prod.world",
   account_authentication = "DATABASE_AUTH",
   username = "username",
   password = "password",
   save_username_password = "SAVE_USERNAME"
)
mxd = arcpy.mapping.MapDocument(mxd_file)
arcpy.mapping.ExportToPDF(mxd, pdf_file)

I am unsure of what I should use for the parameters server and service? When I look at the connection properties of an SDE file created from ArcCatalog it doesn't have those options. Instead it has the single option "Instance" with the value "sde:oracle11g:prod.world", and "Oracle" selected in the "Database platform" dropdown instead of "Application Server".

Attempt 3: Somehow make the MXD remember the username and password. Unsure of how to do that, though.


texts: Copying Domain and Subtype from File Geodatabase?

texts:
Is there a way to copy the domain from one Dataset to another?

I have a series of layers in Dataset A and I want to split the layers in to tiles using the Intersect tool; saving them in Dataset B.  When I do this, all the domain information is lost.


texts: Vector tools missing in QGIS 1.8 (Ubuntu)

texts:
I installed QGIS 1.8 on an Ubuntu 12.04 Precise Pangolin virtual machine.
everything seems to work fine, but I don't get any of the vector or raster tools that (I assume) are part of the standard installation (e.g. analysis, research, geoprocessing tools)
anyone ran into same problem?
Is there a way to get these without uninstalling and re-installing the program?


texts:
I'm using Raster>Extraction>Contours... I specify the input Raster, use the Select dialog to type a name and pick the DXF file type from the pulldown. But instead of a DXF, I'm getting a folder named [NAME].dxf with a shapefile inside that folder. I remember having this problem before and got around it somehow.

Am I doing something wrong, or what's going on?


texts: What does "Computed -srcwin falls outside raster size" mean when clipping?

texts:
What does the following mean, when I try to clip and save the result as a sqlite or other format?

  Computed -srcwin falls outside raster size of 7807x4771


texts: Exporting point pattern fill as vector data in PDF

texts:
In QGIS 3.2, I have a map with parks/other green areas to be displayed as a simple vector point pattern fill. No matter what I do, the pattern fill is always exported to PDF as raster data, which is pretty annoying.

Any suggestions how to change this behaviour?


texts: How to georeference a map in UTM / WGS 84 (geographic) to WGS84 / UTM?

texts:
I have scanned UTM WGS 84 toposheets, with a geographic grid. I want to georeference them in WGS 84 / UTM (therefore with UTM coordinates). Should I:

Calculate the UTM coordinates from the geographic coordinates and enter them for the GCPs, or
Enter the geographic coordinates for the GCPs and then reproject in UTM?

I think "1" is cumbersome but more logical as there is no back and forth change from planar to geographic projection.


texts:
I have a trips table and a bus_stops table. The trips table has a from_stop_id and a to_stop_id. The bus_stops table has lat and lng information for each bus stop.

How can I order the trips in descending order of popularity, where popularity is defined as the number of other trips in the database that start and end within a particular distance of the given trip? For the moment, I am happy to work with 'as the crow flies' distance between two points.

Neither of the tables is spatial but I am using Postgres and adding PostGIS extension is also possible, so if that's the best way forward then I'd love to hear it.

thanks a million


texts: libtiff.net how to read elevation data from non-tiled GeoTIFF files using C#

texts:
I am using libtiff.net to read elevation data from GeoTiff files and it works fine with tiled files but I can't obtain the elevation data from either stripped or scanline-based files especially the ones downloaded from Nasa.
the code I have tried to read to read elevation data from tiled GeoTiff files and its working so far below
but in stripped or line-based files, tiles are always null
                    height = tiff.GetField(TiffTag.IMAGELENGTH)[0].ToInt();
                    width = tiff.GetField(TiffTag.IMAGEWIDTH)[0].ToInt();

                    modelPixelScaleTag = tiff.GetField((TiffTag)33550);
                    modelTiepointTag = tiff.GetField((TiffTag)33922);
                    modelPixelScale = modelPixelScaleTag[1].GetBytes();
                    modelTransformation = modelTiepointTag[1].GetBytes();

                        tilewtag = tiff.GetField(TiffTag.TILEWIDTH);
                        tilehtag = tiff.GetField(TiffTag.TILELENGTH);
                        tilew = tilewtag[0].ToInt();
                        tileh = tilehtag[0].ToInt();
                        tileSize = tiff.TileSize();

                        for (int i = 0; i < width; i += tilew)
                        {
                            for (int k = 0; k < height; k += tileh)
                            {
                                byte[] buffer = new byte[tileSize];

                                tiff.ReadTile(buffer, 0, i, k, 0, 0);

                                for (int itw = 0; itw < tilew; itw++)
                                {
                                    int iwhm = k + itw;
                                    if (iwhm > width - 1)
                                    {
                                        break;
                                    }
                                    for (int ith = 0; ith < tileh; ith++)
                                    {
                                        int iyhm = i + ith;
                                        if (iyhm > height - 1)
                                        {
                                            break;
                                        }

                                        heightMap[iwhm, iyhm] =
                                          BitConverter.ToSingle(buffer, (itw * tileh + ith) * 4);
                                    }
                                }
                            }
                        }

any advise or idea how to do that?


texts: Finding closest tree in LiDAR data with respect to electricity line?

texts:
I have a LiDAR point cloud.

As you see I have one electricity line. I want to vectorise it and then I want to find the closest tree point to electricity line. How can I achieve?

I am open to all programmes and algorithms.

I thought to create a circle and create a cylinder along the line but it is useless too...


texts: How to update the attribute table of feature using PyQGIS?

texts:
Following code updates the value of field ("name") in console  and shows the result but actual attribute table does not get update .Can anybody suggest the error in code.

layer=iface.activeLayer()

selected_feature=layer.selectedFeatures()

layer.startEditing()

for feature in selected_feature:

    feature["name"]="Test name"

layer.commitChanges()


texts: Adding values below column chart for symbology in ArcMap?

texts:
I select a column chart for symbology in ArcMap.

How do I put the values of the parameter under the columns?

Here is an example of what I want to make:


texts: Placing controls outside map container with Leaflet?

texts:
Can someone tell me how I can place controls outside of the map content with Leaflet?

In this case, I just want to place the layer switch control outside the map object.


texts: Update point coordinates in QGIS after manually moving locations in edit mode

texts: I have imported point coordinates from a database, but I need to manually adjust a lot of the locations.  Is there a way for QGIS to automatically update the point locations in the attribute table after you move your points?


texts: Tiles are missing when request from Openlayer2

texts:
I have configured a Tile server and i followed the Link for OSM Tile server from switch2osm

everything is fine, except when i tried to get the tiles using Openlayers its not serving.

But if i give the same url like

http://192.168.1.127/osm_tiles/3/0/4.png

i can see the tiles

but map which is loaded from the Openlayes is not showing the tiles

when i check from the developer console of the chrome the response is HTTP code 200.

Anything i need to check?


texts: Update a field for similar values in the second field based on third field value

texts:
So I have the table like this:

City Distance NearFeature Comments
A    90.5     131025
A    0        781446
B    15.5     131057
B    42.2     140521
C    21.7     784116
C    0        131025

Distance field contains distance value from the City to the NearFeature. Zero in the Distance field means that the City is completely inside the NearFeature.

What I want to accomplish is to populate the Comments field with the value from the NearFeature for Cities which are inside any NearFeature…. Like this:

City Distance NearFeature Comments
A    90.5     131025      City within 781446
A    0        781446      City within 781446
B    15.5     131057
B    42.2     140521
C    21.7     784116      City within 131025
C    0        131025      City within 131025

For B the Comments field is empty because it is outside any NearFeature.

Here is what I’ve done so far. This code updates the table, but only for rows which have 0 in the Distance field. What do I need to change to make it work in the way I need?

cities = arcpy.da.SearchCursor(Path+"\\city.shp", "City")
citynames = []

with cities as rows:
    for row in rows:
        name = row[0]
        if name not in citynames:
            citynames.append(name)
del row, rows

fields = ["City", "Distance", "NearFeature", "Comments"]
with arcpy.da.UpdateCursor(Path+"\\city.shp", fields) as rows:
    for row in rows:
        for i in citynames:
            if str(row[0]) in i and row[1] == 0:
                row[3] = 'City within '+ row[2]
                rows.updateRow(row)
del row, rows


texts:
I'm new to GIS, in PostGIS, and I'm struggling with the SRID concept.

How do you choose a SRID for a database column? This is making me crazy. Why do I need to choose a SRID? What's its meaning?

Picking a SRID is it the same for geometry and geography? Setting a SRID for a geometry field seems senseless, because it should be just on a 2D plane.

When I find lat/lon on the internet or geocoding services, which is the SRID?

Doesn't a lat/lon coordinate make sense without a SRID?


texts: Google Maps has inconsistent zoom level / scale ratio

texts:
At the same zoom level, the scale bar of a map from Google Maps API (v3) can have different sizes.
Here are two screenshots I took at the same zoom level on the same map at two different locations.

Indianapolis

Cambridge (UK)

You can try it yourself on this page (go to the second map, the one with the scale bar).
http://tech.reumer.net/Google-Maps/Demo-Google-Maps-V3/show-all.html

In my application, I have two maps and I want these maps to always have the same scale. I started out by forcing them to be at the same zoom level (with some js listeners), but it looks like it's not going to be enough.

Is there a way to force same scale maps or a specific scale (as opposed to a specific zoom level) on Google Maps API? Which competitors offer this feature?


texts: Error while reprojecting CSV-file with pyproj

texts:
I am trying to reproject a csv-file and here's the error i get.

  File "C:\Users\jonat\Anaconda3\lib\csv.py", line 98, in fieldnames
        self._fieldnames = next(self.reader)

  Error: iterator should return strings, not bytes (did you open the
  file in text mode?)

import csv
import pyproj
from functools import partial
from os import listdir, path

#Define some constants at the top
#Obviously this could be rewritten as a class with these as parameters

lon = 'L' #name of longitude field in original files
lat = 'B' #name of latitude field in original files
f_x = 'x' #name of new x value field in new projected files
f_y = 'y' #name of new y value field in new projected files
in_path = r'C:\H4G_Geomatique\Client\2020\Longueuil\JulienLord\\' #input directory
out_path = r'C:\H4G_Geomatique\Client\2020\Longueuil\JulienLord\\' #output directory
input_projection = 'epsg:4326' #WGS84
output_projecton = 'epsg:2950' #MTM8N

#Get CSVs to reproject from input path
files= [f for f in listdir(in_path) if f.endswith('.csv')]

#Define partial function for use later when reprojecting
project = partial(
    pyproj.transform,
    pyproj.Proj(init=input_projection),
    pyproj.Proj(init=output_projecton))

for csvfile in files:
    #open a writer, appending '_project' onto the base name
    with open(path.join(out_path, csvfile.replace('.csv','_project.csv')), 'wb') as w:
        #open the reader
        with open(path.join( in_path, csvfile), 'rb') as r:
            reader = csv.DictReader(r)
            #Create new fieldnames list from reader
            # replacing lon and lat fields with x and y fields
            fn = [x for x in reader.fieldnames]
            fn[fn.index(lon)] = f_x
            fn[fn.index(lat)] = f_y
            writer = csv.DictWriter(w, fieldnames=fn)
            #Write the output
            writer.writeheader()
            for row in reader:
                x,y = (float(row[lon]), float(row[lat]))
                try:
                    #Add x,y keys and remove lon, lat keys
                    row[f_x], row[f_y] = project(x, y)
                    row.pop(lon, None)
                    row.pop(lat, None)
                    writer.writerow(row)
                except Exception as e:
                    #If coordinates are out of bounds, skip row and print the error
                    print (e)


texts: Using 'gdalwarp' in Python script for Grass

texts:
I'm writing a Python script that aims to re-project various LANDSAT tiles into a common projection and then import the tiles into Grass. The script at present is shown below:

#!/usr/bin/python

import os
import sys
import glob
import grass.script as grass

def import_tifs(dirpath):
    for dirpath, dirname, filenames in os.walk(dirpath):
        for tile in filenames:
        if tile.upper().endswith('.TIF'):
                full_path = os.path.join(dirpath, tile)
                name = os.path.splitext(tile)[0]
        rename = name[:-4] + '.' + name[-2]
        new_path = os.path.join(dirpath, rename)
        re_project = 'EPSG:32643'
                grass.message('Importing %s -> %s@%s...' % (tile, name, dirpath))
        grass.run_command('gdalwarp',
                  '-t_srs',
                  re_project,
                  full_path,
                  new_path)
                grass.run_command('r.in.gdal',
                                  flags = 'e',
                                  input = full_path,
                                  output = name,
                                  quiet = True,
                                  overwrite = True)

def main():
    if len(sys.argv) == 1:
        for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
            import_tifs(directory)
    else:
        import_tifs(sys.argv[1])
if __name__ == "__main__":
   main()

I'm having real trouble using the 'gdalwarp' command to reproject the LANDSAT tiles. I keep getting the following error:

Traceback (most recent call last):
  File "C:/Users/Simon/Documents/import_landsat3.py", line
76, in <module>
    main()
  File "C:/Users/Simon/Documents/import_landsat3.py", line
73, in main
    import_tifs(sys.argv[1])
  File "C:/Users/Simon/Documents/import_landsat3.py", line
31, in import_tifs
    new_path)
  File "C:\GRASS64\etc\python\grass\script\core.py", line
189, in run_command
    ps = start_command(*args, **kwargs)
  File "C:\GRASS64\etc\python\grass\script\core.py", line
167, in start_command
    args = make_command(prog, flags, overwrite, quiet,
verbose, **options)
  File "C:\GRASS64\etc\python\grass\script\core.py", line
124, in make_command
    raise ScriptError("'-' is not a valid flag")
grass.script.core.ScriptError: "'-' is not a valid flag"

I'm assuming this is related to '-t_srs' but I haven't been able to work out a way around the problem - can anyone suggest how to do this?

Many thanks.


texts: What are the best practices for importing and aggregating tracks (in my case in 'plt' format) in a PostGIS table, grouped by user. Are there any tools/db models/tutorials available?


texts:
I have this layer from Geoserver which I overlay in PNG Format on my leaflet map:

/////  Geoserver Layers in WMS format
var pop_tot = L.tileLayer.wms("http://localhost:8081/geoserver/cite/wms", {
    layers: 'cite:vnm_polbnda_adm3_2014_pdc',
    format: 'image/png',
    transparent: true,
    version: '1.1.0',
    attribution: "test"
});

I do have also this part of code to display the labels when I click on my polygons.

///// Add base layers + layers geoserver

    L.control.layers(baseLayers,overlays).addTo(map);

            /// Popup limite_inter_district
                var owsrootUrl = 'http://localhost:8081/geoserver/cite/wms';
                var defaultParameters = {
                    service : 'WFS',
                    version : '2.0',
                    request : 'GetFeature',
                    transparent: false,
                    typeName : 'cite:vnm_polbnda_adm3_2014_pdc',
                    outputFormat : 'json',
                    format_options : 'callback:getJson',
                    SrsName : 'EPSG:4326'
                };
                var parameters = L.Util.extend(defaultParameters);
                var URL = owsrootUrl + L.Util.getParamString(parameters);
                var ajax = $.ajax({
                    url : URL,
                    dataType : 'json',
                    jsonpCallback : 'getJson',
                    success : function (response) {
                    L.geoJson(response, {
                      style: function(feature) {
                        return {stroke: false, fillOpacity: 0.0};
                      },
                            onEachFeature: function (feature, url) {
                                popupOptions = {maxWidth: 250};
                                url.bindPopup("Adm3 Name: " + feature.properties.adm3_name
                                    + "
Total Population: " + feature.properties.pop
                                    ,popupOptions);
                            }
                        }).addTo(map);
                    }
                });

I would like the Pop Up to be active only when my layer pop_tot is selected. How I can do that?


texts:
I'm working with a large dataset of sentinel2 images (jp2 format) and I want to convert them to Gtif format using the comand gdal_translate. I want to create a loop to perform this task. I'm working on Terminal Mac OS 10.12.3 and my GDAL version is 2.1.2, released 2016/10/24.  When doing it one by one, it works using this code

gdal_translate -of GTiff /path/to/input/S2_b1.jp2
/path/to/output/S2_b1converted.tif

By the moment I have tried:

cd /path/to/working/directory for f in *.jp2 do gdal_translate -of GTiff

Any idea?


texts: Is there a way to customize the Options dialog in QGIS 2.14.6?

texts:
I'm creating a deployment package to update QGIS from V2.8 (wien) to V2.14 (essen (LTR)).

In QGIS 2.14 there is a new "Advanced" section under "Settings > Options", see screenshot:

Can this menu item be removed from the Options menu?  I don't want to completely disable the "Options" menu as it is used, however I don't want the users access this advanced set of options.


texts: How to edit data from a feature in Heron

texts:
I'm new to Heron and I would like to know how to add a feature to a layer and to create a custom form to edit data from a feature.

I've tried a few examples, I can select a feature but can't change the data.


texts: Seeking MOOCs (online training courses) for remote sensing?

texts:
I'm looking for a MOOC (Massive Open Online Course) for remote sensing.

Does anyone know if there is currently running or if is there  any change to begin in the future any MOOC for remote sensing?


texts: Performing Multicriteria Analysis Using Distance Along Road Network?

texts:
There is the possibility that I will be working on a project which will look at thousands of properties in a region and measure their accessibility to different services (schools, parks, hospitals, etc).  The idea is that this will highlight properties which are best serviced and properties which need better services.  I've decided that a raster-based multi-criteria analysis would be ideal for this task.

For previous MCA which I have performed, I've always used the Euclidean Distance for calculating distance from objects (water, forests, etc).  However, for this example, this doesn't really make a lot of sense since the services would most likely be accessed via the existing road network.

In the past, I have also performed an accessibility analysis using the Two Step Floating Catchment Area Method (here and here) to determine accessibility to doctors.  I used the road network (with Network Analyst Extension) for this process.  The accessibility was applied to a suburb level, not to an actual property. I was thinking I may need to create an attribute for accessibility to all of these services for each property.  I could potentially weight them as another attribute.  However, an answer to Performing Multicriteria Analysis using QGIS? states that a raster approach is more ideal.

Is there a way to incorporate the accessibility to these services in the form of a raster?

Note, we have not yet got the job and I do not have any data to post.  Let me know if more clarification is needed.

I have ArcGIS and could potentially use the Spatial Analyst and Network Analyst Extensions.


texts: How does ArcGIS project-on-the-fly in an Export Data task, to a different datum and custom projection?

texts:
I need all my data in an equal area projection to for an analysis that will calculate area at the North American Scale. The projection I want is custom, it has the following settings. I will be using the Auxiliary Sphere setting  2 or 3 to improve the accuracy of the area calculations.

Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area
Projection: Lambert_Azimuthal_Equal_Area
False_Easting:  0.00000000
False_Northing: 0.00000000
Central_Meridian:   -100.00000000
Latitude_Of_Origin: 45.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_Sphere_ARC_INFO
Datum:  D_Sphere_ARC_INFO
Prime Meridian:     Greenwich
Angular Unit:   Degree

The data I need to project and transform is in a couple of different formats:
(1) ITRF92 Mexico
and
(2) Canada_Albers_Equal_Area_Conic with GCS_North_American_1983
and
(3) USA_Contiguous_Albers_Equal_Area_Conic_USGS_version with GCS_North_American_1983

I have tried the project tool but there is no valid transformation in ESRI's transformation list that is valid to transform any of the three data sets into my custom projection and so it fails. However if I export the data that needs to be transformed and select the data frame which is set to my custom projection, it "seems" to work. How? Is project-on-the-fly through the Export Data operation reliable? Accurate?


texts: How to tilt/change angle of view/camera from ortho in MapServer or OpenLayers?

texts:
I am wondering if there is a way of changing the angle of the camera in MapServer from Ortho to a 60 degree (+-10) angle.

Perhaps this is an extensive modulation of Google Maps?

I could imagine that with some transformation of the map, this could be achieved by changing, not the angle of view, but the angle of map. So instead of looking at a flat 2d map, then the angle of map could be made by changing the "latitude" of the spheral globe. I then assume distortions of the projection would be an issue, even on a local scale.


texts: Zoom to extent of ajax-loaded vector layer in OpenLayers 3

texts:
I'm loading GeoJSON-data as result from a query to a PostGis-database via AJAX. The feature show up on the map, but I cannot get my map to center on the loaded feature.

var myLayer = new ol.layer.Vector ... // from returned GeoJSON

searchResults.setLayers(new ol.Collection([myLayer]));  // add the layer to a layergroup which is already bound to the map
myLayer.setVisible(true);  // the layer shows up keeping the initial mapview
var extent=parzellen.getSource().getExtent();
alert (extent);  // at this point alert shows me Infinity,Infinity,-Infinity,-Infinity
map.getView().fit(extent, map.getSize());  // this never happens

As shown in the code, propbably due to the fact, that myLayer is called via AJAX the code is executed before the source is completely loaded. Therefore alert(extent) returns Infinity,...

I found, that I have to wait until the source is ready. But in my case this isn't working and the layer doesn't show anymore on the map?

   myLayer.on('change',function(e){
      if(myLayer.getSource().getState() === 'ready') {
        if(myLayer.getSource().getFeatures().length>0) {
          searchResults.setLayers(new ol.Collection([myLayer]));
          myLayer.setVisible(true);
          map.getView().fit(myLayer.getSource().getExtent(), map.getSize());
        }
      }
    });


texts:
I am a recent convert to QGIS.

I have surveyed (and corrected) DGPS (vector point) data collected to OSGB36 (British National Grid); these positions match acceptably well with the OSGB36 map layers to which I have access when imported as CSV and converted to SHP in QGIS, projected using EPSG:27700 (~1m error, although error could be in map layer which is suspect).

I'm attempting a four-stage export from QGISv1.7.4 to Google SketchUp (SU) via Google Earth (GE), as follows:

stage 1: In QGIS, Save as../KML with CRS set to WGS84 EPSG:4326 (or ETRS89);
stage 2: Open KML in GE;
stage 3: Save as KMZ;
stage 4: Open in SU.

Unfortunately, the exported WGS84 positions in the KML file are in error by ~8m - I've checked them against OSTN02 conversion using tool on OS website , I've also checked GE projection (by eye: converting GE lat/long to OSGB36 and comparing with OSGB36 vector map layer).

SU fails to open the KMZ file as well, but that's not a QGIS issue, I imagine.

So the evidence seems to point towards an error in the transformation from OSGB36 to WGS84 in stage 1, i.e. with QGISv1.7.4.

Do you have any suggestions as to where the fault in my process might lie, or to an alternative approach?


texts: Cut polygon results do not equal to the source data

texts:
I have a polygon that I would like to cut as two parts. After using the cut polygon tools, it has divided into two parts. However the total area of that two parts does not equal to the original polygon.

Example: Total area of two parts is 239.535215 m^2, Area of original polygon is 239.535226 m^2

How can I get the cut polygon result which is exactly equal to the original polygon?


texts: reading rasters in r using terra package

texts:
Using the R terra package for the first time after many years using the raster package. Stuck on the basics.
Trying to read DEM raster into r and reproject it. The raster is in ESRI grid format and can be downloaded here.
I previously asked this question on r-sig-geo.
And had working solution with raster
dem_raster<-raster("/Users/mypath/w001001.adf")

Attempting the same with terra::rast
dem_terra<-rast("/Users/mypath/w001001.adf")

And plotting the results side by side...

Printing the details for the two rasters
dem_raster

class      : RasterLayer
dimensions : 8827, 10224, 90247248  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : -5762000, 4462000, -3920000, 4907000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs
source     : w001001.adf
names      : w001001
values     : -76, 5930  (min, max)
attributes :
         ID COUNT
 from:  -76     2
  to : 5930     1

dem_terra

class       : SpatRaster
dimensions  : 8827, 10224, 1  (nrow, ncol, nlyr)
resolution  : 1000, 1000  (x, y)
extent      : -5762000, 4462000, -3920000, 4907000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs
source      : w001001.adf
name        : VALUE
min value   :   -76
max value   :  5930

How do I get rast to read in the ESRI grid values correctly?


texts: Controlling output linewidth linestyle using gdal_rasterize shapefile shp to Geotiff

texts:
I have successfully combined a shapefile (.shp) and a GeoTiff using gdal_rasterize.

The shapefile displays as one pixel wide in the GeoTiff, but I need it to be wider and more visible, say, 5 pixels wide and red.

Can anyone suggest an alternate way to achieve this? This needs to run as a script and not use desktop gis software.


texts:
I've imported a shapefile from the Corine Land Cover dataset in PostGis. I'm working on a software which will allow for diverse simulations on the territory.
This software integrates a function to connect in a database PostGis, my shapefile CLC when imported into PostGIS has no classification of the type of ground in the display.

However, when I import directly to the shapefile "raw" without passing through PostGis, the type of ground is shown well on the software.

This software already integrates the table with the codes CLC but it does not seem to work with tables of PostGis...

Thus it is a problem at the level of PostGis? Nevertheless in the table, I have well the column which defines the type(chap) of land use.


texts:
I have several tiles in a catalog covering bays along the shore. There are some hits in the water but they are sparse. Ideally I think I should use the is_concave = TRUE as an option in grid_terrain(), but it was taking so long I ended up terminating R.
I looked at the resulting terrain rasters that had the "outside the convex hull" warning - they had a small area of extreme negative values due to the NAs in the water area. Is it possible to remove the NAs?
Data may be downloaded from NS elevation explorer (search for 284_4955).
library(lidR)

myLas <- readLAS("284_4955_201901.laz", filter = "-drop_withheld -drop_overlap -drop_z_below -1")
x <- grid_terrain(myLas, 1, tin())
#> Interpolation of 89 points outside the convex hull defined by ground points (outside the triangulation) failed and returned NAs.

x
#> class      : RasterLayer
#> dimensions : 1000, 1000, 1e+06  (nrow, ncol, ncell)
#> resolution : 1, 1  (x, y)
#> extent     : 284000, 285000, 4955000, 4956000  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=20 +ellps=GRS80 +units=m +vunits=m +no_defs
#> source     : memory
#> names      : Z
#> values     : -21474836, 106.02  (min, max)

2 diagonal black lines are negative raster values which I think come from the NAs.


texts: Sentinel 3 OLCI/SLSTR granule footprints shapefile?

texts:
Basically the following question:
Finding Sentinel tile for specific Long/Lat coordinate

However, that question was regarding Sentinel 2 MSI, for which a UTM tiling grid is made available by the ESA.

For a given study area, I would like to find out which Sentinel 3 OLCI (and possibly SLSTR) granules are needed. Is anyone aware of a file describing the footprints of these granules?


texts: How to add a feature on a ArcGIS Feature Service with a Geocortex Workflow

texts:
Do you know the correct Geocortex workflow to add a feature to a layer of a ArcGIS Feature Service?

I'm using the Geocortex Workflow Designer 4.9.2.46. I know it's easier to add a feature with the Geocortex Workflow Designer 5.

The sequence of activities I'm currently using is the following:

Display Capture Geometry
Display Form
Add To Dictionary (here I'm adding the Captured Geometry to the dictionary returned by the form)
Convert the dictionary to a DataTable
Convert the DataTable to a FeatureSet
AddFeatures

The feature is effectively created but the geometry attribute is missing. I don't know how to update this geometry field. In the database, this geometry field is called shape.


texts: Use UpdateCursor to update a field based on text in other field

texts:
I am trying to find the correct syntax for the UpdateCursor where I look for part of a text from another field.  I have many values of PCD zoning codes (PCD 03, PCD 04, PCD 05, etc).  I want to calculate my "ZONING_DESC" field based off these PCD values in my "ZONING" field where all PCD values get calculated to "Planned Commercial Development".  I know that '==' operator is finding exact values.  Is there an operator for finding like values?

Here is my Code:

with arcpy.da.UpdateCursor(Zoning, ["ZONING", "ZONING_DESC"]) as cursor:
    if row[0] == 'PCD%':
        row[1] = 'Planned Commercial Development'
    cursor.updateRow(row)


texts: How to create an index map of atlas features in QGIS 3?

texts:
I realise that this question has been asked before (Making overview with linked detailed maps in QGIS?), however I am still unable to figure out how to create an index map of atlas features in the QGIS map composer.

Basically, I have a line and I have five atlas maps showing the extent of the line, and to show where which map corresponds to the line a I need five overviews on a map with numbers on them, which will look similar to below:

I've created five overviews, but as all of these must be set to correspond to the main map frame, all of them show the overview for which ever atlas feature is selected at the time, rather than appearing as separate index tiles.

Can anyone tell me where I am going wrong?


texts: Create custom locator file - what license level of ArcDesktop do I need?

texts:
I would like to create custom locators in ArcDesktop to circumvent this problem:
US One Address locator style missing city, state and zip code fields?

So I am going to buy a copy of ArcDesktop, but I think there are 3 levels, so I'm not sure which one to buy?  ArcView, ArcEditor or ArcInfo.

Which level allows you to create custom locators, aka .lot files?


texts: Displaying multiple ArcMap documents in one file

texts: I have a project in which I'll have 22 different maps that I'll need to display on a single 3' x 4' poster. Each map is an individual .mxd file in ArcMap. How can I best accomplish putting these 22 maps in a single visual array?


texts: Loaded vector layers from SQL Server 2014 but data is not visible on the map in QGiS 3.8.3 - Zanzibar

texts:
This is odd - I've successfully connected to my spatial sql server 2014 instance, and loaded two point layers and one Polygon layer into QGIS.  I see nothing on the map.  I zoom to the layer, I know it goes to the right location and zoom level.  I use the Info option, click the map, and it doesn't find anything.  Very puzzling.  The Style appears to be fine.  I don't know what else to check.  All looks good configuration wise - I just don't see anything.

I then selected each of the three layers one at a time, and used the Layer drop down menu, SaveAs... and saved the layer as a shape file at my local computer.

When I load the three Shape files that I've just created by saving the SQL Server layers as shape files using QGIS, these shape files display and act perfectly.

Has anyone else experienced this?  I don't know why I cannot successfully work with SQL Server 2014 spatial layers.

For the record, all my layers use the EPSG:2234 CT State Plane CRS, in case that may be relevant.  The geometries in the database do specify this CRS, and QGIS displays this as well, and it is set correctly in the saved shape files.


texts: Recover failed raster calculator operation QGIS

texts:
Users can recover their raster calculator algorithms via the "results" panel in ArcMap (i.e., analogous to viewing one's "history"). Is there a similar function in QGIS?

In other words, is there a way to recover raster calculator equations after they have already been (unsuccessfully) run?


texts: What's your view on the direction of GIS?

texts:
I am a long standing GIS developer and consultant.  I am often asked, is there still a place for Desktop GIS, to which I, currently, answer, of course; if there are jobs specific to individuals, you shouldn't be developing enterprise solutions for them.

However, with the advent of the cloud (for want of a better descriptor), and the advances in delivering Web Services and functionality via the cloud, how long do you think it will be before enough functionality is deliverable via the cloud?  Or do you think, as I do, there will always be a place for desktop GIS?

I am interested, as it is a developing situation and most of my clients are interested in knowing where we (GIS) are heading.  For now, I still believe that the application or resolution of complex GI problems will necessitate the need for Desktop apps and that it will be some time before all the necessary functionality is deliverable via the web AND, I believe it is still not worth the extra development for companies in delivering or developing complex solutions for power GIS users.

What say you?


texts:
I have two polygon layers, and a roads layer that I've made into a network dataset. My goal is to find the closest polygon in one layer to each polygon in the other layer, but instead of just calculating the distance to find the closest I need to use the roads ND to see which of the polygons from the first layer is closest to the polygon from the second layer.
What tool can I use to do this?
I've assumed I can do it with something in the Network Analyst toolset, but I can't figure out what.
From the picture, I want to find the closest green polygon to both of the red polygons using the yellow roads.
I'm still relatively new to this, in my first course on the subject.

Edit: I received this solution on GeoNet, the Esri forum:

Convert the vertices of your polygon features to points (Feature Vertices To Points-Data Management toolbox).
Do this twice: once for the green polygons and once for the red polygons
Calculate an OD cost matrix with the Network Analyst (OD cost matrix analysis layer-ArcGIS Pro).
Load the vertices from the red polygons as origins and the vertices from the green polygons as destinations (or vice-versa).
Use the OD cost matrix to find out the shortest path between any vertex of a green polygon and any vertex of a red polygon.
Use these vertices geographic positions to determine to which green or red polygon the vertices belong.


texts:
I'm having a hard time dealing with geotiff images on a geoserver: I can't seem to get rid of the black 'no data' background around my geotiff image (a map area with a transparent background).

Further clarification: the geoTiff is in a geoserver layer group, and not a layer of its own, and it has the geoserver 'default' raster style on it (nothing but an opacity setting).

By default, my colour image is a black and white one; I can set the colour channels through the layer style, but that strips the alpha info (as I read elsewhere).

Setting the 'transparent colour' (to #000000) in the geoserver layer specifications doesn't seem to work, either.

I also see a place in the layer-specific page where we can set the channels/bands (there are five in mine, but the first four are rgba, verified by gdalinfo) to certain values (range, 'unit' (?)), but I can find no documentation on this.

Can anyone enlighten me?


texts: How to retrieve the per-feature value rendered when IGeoFeatureLayer.DisplayAnnotation is true (ArcObjects with C#)?

texts:
In another area of my application I'd like to display the per-feature label value rendered when Labeling is turned on in ArcMap, i.e. when IGeoFeatureLayer.DisplayAnnotation becomes true.

I have a list of items of type IFeature from say myFeatLayer.  I am able to retrieve the layer-level labeling properties, including text formatting, with:

IGeoFeatureLayer geoLayer = (IGeoFeatureLayer)myFeatLayer;
IAnnotateLayerPropertiesCollection anno = pGeoFeatureLayer.AnnotationProperties;
IElementCollection iecp;
IElementCollection iecup;
IAnnotateLayerProperties props;
anno.QueryItem(0, out props, out iecp, out iecup);
ILabelEngineLayerProperties leprops = props as ILabelEngineLayerProperties;
ITextSymbol text = leprops.Symbol;

What I don't know is necessarily how to query an IFeature, and with what (leprops.Expression?), to retrieve the actual value of what will be rendered near each feature by the labeling engine.

Thanks in advance.


texts: Copying a field (column) value from a Voronoi polygon into the point inside that polygon

texts:
I'm relatively new in GIS and I am using QGIS. I'm trying to find a way to copy a field value from a voronoi polygon layer into each point that is currently inside the said polygon.
Each voronoi polygon contains a unique population field value. Also, each polygon contains exactly 10 points distributed randomly within. What I intend to do is to copy the population values from each polygon into the points within the same polygon then divide them by 10 - essentially distributing the population to every point.


texts:
I am trying to load a WMS into my C++ program to use it as a GDALDataset. The problem is that I am able to load the WMS but every time I am doing a GetRasterCount() I get 0; and every time I am trying to do a GetRasterBand(0) I get "Erreur de segmentation (core dumped)".
I was doing the same thing for WFS (GetLayerCount() and GetLayer(0)) and it was working fine.
Here is my code:
GDALAllRegister();
GDALDataset *poDS = static_cast<GDALDataset*>(
poDS = (GDALDataset*)GDALOpenEx("https://geowebservices.stanford.edu/geoserver/ows?SERVICE=WMS&request=GetCapabilities", GDAL_OF_RASTER,NULL,NULL,NULL ));

if( poDS == NULL )
{
    std::cout << "Dataset is empty" << std::endl;
    GDALClose( poDS );
    exit( 1 );
}

std::cout << "THE DATASET CONTAINS THIS MANY RASTER BANDS : " << poDS->GetRasterCount() << std::endl;

for (int i = 0; i < poDS->GetRasterCount(); i++){

    GDALRasterBand  *band;
    band = poDS->GetRasterBand(i);
    std::cout << band->GetXSize() << std::endl;
}


texts: Setting all values in column to be same using QGIS?

texts:
Using QGIS, I want to fill a whole column with the same "word".

Is there a faster way or an expression I can use in the field calculator instead of copy/paste?

In some layers I have more than a thousand attributes.


texts: gpx conversion tools in qgis keeps asking for GPSBabel

texts:
I want to perform a simple task: show the waypoints of a GPX "track" as dots, not as a set of connected lines in QGIS ( Conversion -> Waypoints from a track ). I use files, not any particular USB connected device. Just GPX standard files. QGIS keeps asking for a "GPSBabel" module. Did I miss something? The same data is imported as a "track" without problem.

( I'm using the version 1.8 of qgis )


texts: How to Georeference a tiff with correct coordinate system

texts:
I am having difficulty in georeferencing an existing tiff (already has georeference information - but is out of alignment). Using the georeference plug-in I can create the map's GCP in the 4 corners and run the transformation (using WGS84). However, even trying other transformations I get a map that is not georeferenced to the Lat & Long. The Lat / Long of the original map shows 26.00, -25.00 (top left corner), the georeference after transformation is 974, -747.

How can I get the transformation to reflect Lat / Long as per the original georeference as all I am trying to achive is a re-aligned map?


texts: How can I project a raster DEM and not get serious distortion in the resulting DEM?

texts:
Please see this link to see a fine example of what I am dealing with.

Distortion in projected raster

My specific situation is that I have downloaded DEM data from the National Map (USA) and it comes in GCS North American 1983 with no linear unit. Without a linear unit you cannot calculate slope, so in order to calculate the slope I project the data in the appropriate UTM zone (using both ArcGIS 10.4 and ArcGIS 10.2, Project Raster tool) and what I get is a DEM that looks like it was made of clay and someone smashed it up against a chain link fence (just like the examples in the aforementioned link).

I see that the linked question above was closed for not being specific enough, but nobody really provided any kind of answer either, and to me the question seemed pretty clear. Maybe we are both going about this the wrong way, is there a better way to project a raster or assign a linear unit to a GCS dataset?

I would use a different dataset, but there really isn't one and I feel like this data should be usable since it is the National Elevation Dataset.


texts:
I want to use OpenTopoMap as background map in Leaflet.

I can't find information about the colour encoding of heights to create a legend.

Does anyone know the correspondence between colours and height intervals?


texts: How to create foreign keys and cardinality in a Geopackage with QGIS or FME? It's easy to get the tables, but how can I set the relations between the tables?


texts: GIS shapefile of corporate headquarters located in the US?

texts: I am conducting MSA (Metropolitan Statistical Area) level analysis where I am tracking the number of corporate headquarters in the top 50 MSAs from 2000 to 2014. I am looking for a shapefile like this for both 2000 and 2014 (or close to 2000 and 2014). Has anyone used such shapefiles and if so, did you find them online or have to construct them on your own? If there is a shapefile for all corporate headquarters in the US for 2000 and 2014, I can do select by location with a MSA shapefile. Or if I can get a list of corporate headquarters by MSA for 2000 and 2014, I can just join the data file to the MSA shapefile.


texts: Displaying a dynamic KML file on Google Maps as it is on Google Earth

texts:
how can i display my kml files through network link on Google Maps instead of showing it on Google Earth.

I have tried the step mentioned on this link

https://developers.google.com/kml/articles/pagesforkml

but stucked at step 3 where it is saying for uploading kml file as no option is available there.

Actually ,I'm uploading KML maps because , my KML file size is more than 2 MB and as i refresh it on google earth , it got hanged..

So what i need is a network link on Google maps which will show my my.kml file placemarks and this network link is having a refresh interval time ..my.kml file is generating at every 10 seconds by a java program...

my network link code is :
   // pasted only the network link part

<NetworkLink>
<Link>
<href>/home/placeMarkers.kml</href>
<refreshmode>onChange</refreshmode>
<refreshInterval>60</refreshInterval>
<viewRefreshMode>onStop</viewRefreshMode>
<viewRefreshTime>20</viewRefreshTime>
 </Link>
</NetworkLink>
</kml>


texts: How best to use shapely to write a dict to a shapefile

texts:
I have created a shapely object (i.e. the_pt) for my Points & Multipoints to do a buffer analysis against another shapely object (i.e. city_pt) using:

  return_shapes = []
  for pt in points:
            # create a shapely obj of all the points
            the_pt = shape(pt['geometry'])
            # what type of geometry is it?
            geom = pt['geometry']
            # check if it is a multipoint or point
            if geom['type'] == 'MultiPoint':
                # i.e. <class 'shapely.geometry.multipoint.MultiPoint'>
                pts = []
                for m_pt in the_pt:
                    city_pt_buffer = city_pt.buffer(dist)
                    # city_pt is a shapely obj
                    if m_pt.within(city_pt_buffer):
                        pts.append(m_pt)
                return_shapes.append(MultiPoint(pts))
            # check if it is a point geom
            elif geom['type'] == 'Point':
                    city_pt_buffer = city_pt.buffer(dist)
                    if the_pt.within(city_pt_buffer):
                        return_shapes.append(the_pt)
        return return_shapes

If I evaluate pt - I have a dict type with the geometry and properties needed to map to shapefile...

 <type 'dict'>
    {'geometry': {'type': 'Point', 'coordinates': (-13649075.348882342, 6193380.937753927)}, 'type': 'Feature', 'id': '7', 'properties': OrderedDict([(u'OBJECTID', 68547), (u'ISO', u'USA'), (u'SITE_NAME', u'Anacortes')])}

however what I really want is 'return_shapes' to be written to shapefile (with all the properties...Is there anyway to build return_shapes (a list) into something I can write (map) into a shapefile format?


texts: Covert /z/x/y.jpg map tiles to something that can be provided using the WMS specification

texts:
We have generated some maps in our own software solution where the output is a map stored in XYZ tiling format.

/zoom/x/y.jpg

I got a request if I could provide it in a WMS.  I am wondering if anyone has some experience on how to in shortest amount of time can convert and setup a WMS that can service these images.

The reason is that the map client dont support XYZ and is bound to use normal WMS only.

We would prefer a windows solution but anyhting will do.


texts: QGIS: Python script to batch clip vector layers and set their CRS

texts:
Here below the sript I'm trying to run.
clip_zone is a polygon and in my input_folder, I have a point layer and polyline layer.

#Tell Python where you will get processing from
sys.path.append('C:\Program Files (x86)\QGIS Brighton\apps\qgis\python\plugins\processing')

#Reference the algorithm you want to run
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

#Define the workspace
input_folder = "C:\Users\id846890\Documents\FAN\Design\Test_clip\input\ "
output_folder = "C:\Users\id846890\Documents\FAN\Design\Test_clip\output\ "
clip_zone = QgsVectorLayer('C:\Users\id846890\Documents\FAN\Design\Test_clip\FH.shp','FH', 'ogr')
crs=clip_zone.crs()
crs.createFromId(31370)
clip_zone.setCrs(crs)

#Clip the shapefiles
for vectorLayer in input_folder:
    layer_clipped = os.path.join(output_folder, vectorLayer.name())
    processing.runalg('qgis:clip', vectorLayer, clip_zone, layer_clipped)


texts: I'd like to know how GPS can determine if a street is one-way. Does it detect one-way signs or parked car positions if any?


texts: Problem to create composer with OpenLayers plugin in QGIS

texts:
Creating layout in composer I am facing a problem of dislocation of a shapefile and a Google Earth image (which I am using from "OpenLayers plugin") Version 2.4.0 Chugiak (Mac).

In the first image shows the situation of map composer QGIS. Everything is correct and I just wanted the same output, but when I save it as PDF or JPEG or something, the shapefile layer is seen as a bit moved..!

Just want to get the exact same image which I am getting in the composer!


texts: How can I develop a vehicle tracking application on a tight budget?

texts:
I need to develop and GIS application in .net with following Requirement,
Can you please suggest the options how can i start with.

I am new to GIS World, but I have a good knowledge.
I neither developed any GIS application  before!

1) We need to track vehicles on the map, the data recieved from Multiple Controllers connected Over Ethernet.

2) Display of Movements of Multiple Vehicles to be Provided for simultaneous operation.

3) Displaying Latitude, Longtitude, Altiitude, Roll, yaw, Pitch, Velocity with Every Vechicle.

4) Unique Icons for Multiple vehichles (icon should be configurable)

5) Provision for the modes to facilitate selection of vehicles for simultaneous display on Map.

6) Automatically move the Vehicle icon  representing the respective vehicle, to the correct geographic location on map when vehicle change the position.

7) Automatic PAN movement of the area concerned on the digital map during Online Vehicle Tracking.

8) Calculate Actual disptance traveled by vehicle along with Aerial Distance,

9) Calculate relative distance between two vehicles.

10) Digital Map retrieval and display as per the area of operation.

11) Provision for additional icons creation on the map. (These icons should be displayed for futher use)

12) Provision of Information Layer (Peliminary/Detailed information about Terrain etc), 3D ANALYSIS

13) Geo-Referencing Facility. It should be able to digitize it by georeferencing.

14) Import digital masp irrespective of their original/native format.

15) OS Platform independency.

It should support Military grid coodrinate Data.

17) Provision for Defininf Waypoints with positioning Information

18) Multiple Access to various Digital map,image map like geotiff,.jpeg,.dted image

19) Integration of scanned maps, satellite images .

20) Zoom in Zoom out, Zoom in Max Scal1 1:1 and zoom out 1:50000

21) User Specified scale and Dynamic display of scale in the Mapview

22) Exchange of GIS Data for generation of Reports.

23) Import/Export Maps using .dgn (Microstation DGN) format.


texts: Best way (performance wise) to add Attribute fields in PyQGIS on large files

texts:
I have a question regarding adding attribute fields on an existing vector file in PyQGIS. I’m writing a processing script for QGIS, and so far, I use the QGIS field calculator to add new attribute fields and then fill every feature with a value. The script looks something like this:

layer = QgsVectorLayer(myVectorFile,"", "ogr")
Field_1=processing.runalg('qgis:fieldcalculator',layer,'Field_1',1,1.0,1.0,True,Value_1,None)
Field_2=processing.runalg('qgis:fieldcalculator',Field_1['OUTPUT_LAYER'],'Field_2',1,1.0,1.0,True,Value_2,None)
Field_3=processing.runalg('qgis:fieldcalculator',Field_2['OUTPUT_LAYER'],'Field_3',1,1.0,1.0,True,Value_3,None)

Etc.

This works as expected, but it takes a very long time when I run the script, since it runs on very large vector files with a lot of features. I’m looking for a better/more performant way to add fields and file them with a value in PyQGIS. I imagine that maybe it would be possible to only edit the .dbf file? I’m open to suggestions.


texts: Error when ingesting files in Google Earth Engine

texts:
I am trying to ingest a simple point shapefile in GEE from the asset tab of the Code Editor
The program sent me back this error

"Task details: Asset ingestion: users/giuliocst/JesrAdbachPOINT

State: Failed

Started: 1m ago (2018-11-07 16:34:10 +0100)

Runtime: 7s

Id: LVKE25EJBW5AWJ6AVYPDCVKV

Error: File access denied: 'E1yLz1m5tKhVI1Io/Jesr Adbach POINT.shp'."

What I am missing? I am using the same procedure as always


texts: Removing deadlock on PostgreSQL layer by Cadcorp SIS?

texts:
I have just created a new table which I will want available for multiple connections at once. I have never had to consider this before as all tables in the past have been fine for this.

However this new tables seems to be creating a deadlock when two (or more) connections are attempted from Cadcorp SIS 9.

The first connection will open the table fine. The second will get stuck and throws back the following:

2019-04-29 09:57:22.610 BST [7276] ERROR:  deadlock detected
2019-04-29 09:57:22.610 BST [7276] DETAIL:  Process 7276 waits for AccessExclusiveLock on relation 6798789 of database 774436; blocked by process 6168.
                Process 6168 waits for AccessExclusiveLock on relation 6798789 of database 774436; blocked by process 7276.
                Process 7276: select ST_EstimatedExtent('public', 'mastermap', 'geometry') LIMIT 0
                Process 6168: DECLARE C8 NO SCROLL CURSOR FOR select ST_EstimatedExtent('public', 'mastermap', 'geometry')
2019-04-29 09:57:22.610 BST [7276] HINT:  See server log for query details.
2019-04-29 09:57:22.610 BST [7276] STATEMENT:  select ST_EstimatedExtent('public', 'mastermap', 'geometry') LIMIT 0

If I kill the blocking session, it actually allows the second one connection through without affecting the first one (so at this point both work fine). However if you then close one of the Cadcorp sessions and reopen, it just locks again.

How do I remove this type of lock permanently? It is not a dataset that anyone will ever edit, it's only for viewing/tracing purposes and i don't want to have to kill the session every time multiple connections open up.

I have read the docs on locking and unfortunately cannot make sense of them  (I am a SQL novice).

UPDATE AFTER RUNNING SUGGESTION FROM JGH:

INFO:  vacuuming "public.mastermap"
INFO:  index "mastermap_geometry_sidx" now contains 483426581 row versions in 3945914 pages
DETAIL:  0 index row versions were removed.
0 index pages have been deleted, 0 are currently reusable.
CPU: user: 8.43 s, system: 38.37 s, elapsed: 260.44 s.
INFO:  "mastermap": found 0 removable, 4 nonremovable row versions in 1 out of 22441174 pages
DETAIL:  0 dead row versions cannot be removed yet, oldest xmin: 59942
There were 0 unused item pointers.
Skipped 0 pages due to buffer pins, 0 frozen pages.
0 pages are entirely empty.
CPU: user: 9.09 s, system: 38.46 s, elapsed: 261.42 s.
INFO:  vacuuming "pg_toast.pg_toast_776026"
INFO:  "pg_toast_776026": found 0 removable, 9 nonremovable row versions in 1 out of 205740 pages
DETAIL:  0 dead row versions cannot be removed yet, oldest xmin: 59943
There were 0 unused item pointers.
Skipped 0 pages due to buffer pins, 0 frozen pages.
0 pages are entirely empty.
CPU: user: 0.01 s, system: 0.00 s, elapsed: 0.22 s.
INFO:  analyzing "public.mastermap"
INFO:  "mastermap": scanned 30000 of 22441174 pages, containing 644561 live rows and 0 dead rows; 30000 rows in sample, 482156852 estimated total rows
VACUUM

Query returned successfully in 8 min 50 secs.


texts: How to use part of the attribute value as a label in ArcGIS Desktop?

texts:
I have a data (school.shp) about schools. I have attribute values like:

TABIKA KG SENTOSA, TADIKA BANDAR SEMARIANG, KEMAS KG TABUAN, TABIKA KEMAS SANTUBONG, TADIKA KEMAS SERIAN

I need to label it like below::

TABIKA KG SENTOSA - TABIKA

TADIKA BANDAR SEMARIAN - TADIKA

KEMAS KG TABUAN - KEMAS

TABIKA KEMAS SANTUBONG - TABIKA KEMAS

TADIKA KEMAS SERIAN - TADIKA KEMAS

I just need to label TADIKA, TABIKA, KEMAS, TADIKA KEMAS, TABIKA KEMAS, no need to label the full name. Is there any VB or SQL query that can help me to settle this problem?


texts: Testing uppercase vs lowercase field values on the fly using ArcPy & Python?

texts:
I am trying to use the code below to allow testing of one field which contains mixed case against another field which contains uppercase only.
I have read many many posts but they generally tend to look at converting or using Field Calc and I do not want the values to be converted. (only on the fly)
I am using ArcGIS 10.2 Win 7 x64 and Python 27.

    # Refine the selection set to only 'LGA <> LGA_NAME'
    fn99 = "LGA"
    fn98 = "LGA_NAME"
    EXPR1 = fn99.upper() +"<>"+ fn98.upper()
    DEF1 = arcpy.SelectLayerByAttribute_management(fc_lyr2, "NEW_SELECTION", EXPR1)
    print "\n" "Using SelectLayerByAttribute to create a 'NEW_SELECTION' from " + fc_lyr2
    print DEF1.getMessages()

Note I dont get any errors, my code just doesnt test successfully so all values are retyurned and I know there are at least 50 incorrect values from manual checking.


texts: QGIS Renders Buffer Polygons from ArcGIS Feature Class Differently

texts:
I have a polygon feature class stored in an ArcGIS geodatabase. All of the shapes are buffers of a point layer created using the editor menu's 'Buffer' command. In ArcMap 10.3.1 the data looks like this:

Using the ESRI FileGDB driver, I can bring the data into QGIS (tried on 2.8.6 and 2.8.14). However, the drawing of the buffers is completely different:

I understand there may be implementation differences that cause this and I probably can't do much about it. What I'm really interested in is if someone knows what accounts for the difference in rendering of circular polygons between ArcGIS and QGIS, or could point me somewhere (ogr documentation? Qgis documentation?) where I might figure this out.


texts: FileGDB In GDAL Docker container with dot net core

texts:
Currently, I am using this docker image pepperlk/gdal that has FileGDB built in and working via command line.
Details here:
gdalinfo --version
GDAL 3.3.0beta1-b75f14ab12155c6b6f4544a9894c41f2f9a7d241, released 2021/04/20

ogrinfo --formats
<concat for clarity>
  OpenFileGDB -vector- (rov): ESRI FileGDB
  FileGDB -vector- (rw+): ESRI FileGDB

To reproduce...
Running within the container:
MaxRev.Gdal.Core.GdalBase.ConfigureAll();
var drcount = OSGeo.OGR.Ogr.GetDriverCount();
for (int i = 0; i < drcount; i++)
{
    var driver = OSGeo.OGR.Ogr.GetDriver(i);
    Console.WriteLine($"Driver Name: {driver.name}");
}

Yields:
Driver Name: ESRIC
Driver Name: PDS4
Driver Name: VICAR
Driver Name: PDF
Driver Name: MBTiles
Driver Name: BAG
Driver Name: EEDA
Driver Name: OGCAPI
Driver Name: ESRI Shapefile
Driver Name: MapInfo File
Driver Name: UK .NTF
Driver Name: LVBAG
Driver Name: OGR_SDTS
Driver Name: S57
Driver Name: DGN
Driver Name: OGR_VRT
Driver Name: REC
Driver Name: Memory
Driver Name: BNA
Driver Name: CSV
Driver Name: NAS
Driver Name: GML
Driver Name: GPX
Driver Name: KML
Driver Name: GeoJSON
Driver Name: GeoJSONSeq
Driver Name: ESRIJSON
Driver Name: TopoJSON
Driver Name: Interlis 1
Driver Name: Interlis 2
Driver Name: OGR_GMT
Driver Name: GPKG
Driver Name: SQLite
Driver Name: WAsP
Driver Name: PostgreSQL
Driver Name: OpenFileGDB
Driver Name: XPlane
Driver Name: DXF
Driver Name: CAD
Driver Name: FlatGeobuf
Driver Name: Geoconcept
Driver Name: GeoRSS
Driver Name: GPSTrackMaker
Driver Name: VFK
Driver Name: PGDUMP
Driver Name: OSM
Driver Name: GPSBabel
Driver Name: SUA
Driver Name: OpenAir
Driver Name: OGR_PDS
Driver Name: WFS
Driver Name: OAPIF
Driver Name: HTF
Driver Name: AeronavFAA
Driver Name: EDIGEO
Driver Name: SVG
Driver Name: CouchDB
Driver Name: Cloudant
Driver Name: Idrisi
Driver Name: ARCGEN
Driver Name: SEGUKOOA
Driver Name: SEGY
Driver Name: ODS
Driver Name: XLSX
Driver Name: Elasticsearch
Driver Name: Carto
Driver Name: AmigoCloud
Driver Name: SXF
Driver Name: Selafin
Driver Name: JML
Driver Name: PLSCENES
Driver Name: CSW
Driver Name: VDV
Driver Name: GMLAS
Driver Name: MVT
Driver Name: NGW
Driver Name: MapML
Driver Name: TIGER
Driver Name: AVCBin
Driver Name: AVCE00
Driver Name: HTTP

In .net I am missing that installed driver what am I missing??
ogr2ogr is working fine as well as ogrinfo via command line but not available in code.


texts: OpenLayers low resolution and/or shifts in QGIS print composer?

texts:
I'm having trouble with the OpenLayers Plugin of Qgis and the map composer:

If I create an OSM-background layer and if I want to export this, the OSM-Layer looks perfectly all right in the normal qgis program window. But in the map composer and after export the layer has shifted relatively to my other shape layers (EPSG:32633 - WGS 84 / UTM zone 33N).

The second thing is that the output resolution of the exported osm-layer is very, very poor.

A really bad way of getting around this whole trouble would be to increase the screen resolution and make a screenshot of the map composition window of qgis. But I don't think this would be very professional. It also would cause a lot of pain :)

I'm using Qgis 1.8.0-Lisboa under Linux. The openlayers plugin is version 0.92.


texts: ArcGIS Server Java script API error code 500 while using proxy.ashx on one particular server

texts:
I am looking for help with using proxy.ashx in Javascript API on one particular server at a client site. We have a GP service that the proxy is used to access (when request is more than 2048 characters).

I have configured proxy and it works fine on one of our server in house. Javascript API application can send request using proxy.ashx.

In another environment at our client site we get Error 500 -Unable to access the web page. But without the proxy we can access the page.
Here are some test case scenarios we tested.

We have the Javascript application on H01 server and GP Service on H02 server. So we send request from H01 to H02. Both the machines are in the same domain.

Following tests works fine

http://H02/rest/services/MyGPService/GPServer
http://H02/iisstart.htm
http://H02/rest/services/MyMapServices/MyMapService
http://H01/MyApp/proxy.ashx?http:www.google.com

Following tests don't

http://H01/MyApp/proxy.ashx?http://H2/rest/services/MyGPService/GPServer
http://H01/MyApp/proxy.ashx?http://H2/iisstart.htm
http://H01/MyApp/proxy.ashx?http://H2/rest/services/MyMapServices/MyMapService

The "H02" server where AGS is installed is managed by their IT department and we don't have much access to it.

In the tests above, I cannot access any URL from "H02" server using proxy, even the root iisstart.htm. They all give Error 500, the page cannot be accessed.

It looks to me like permission or settings on "H02" server that is causing this issue. All the above tests works fine on our internal server.


texts: GDAL/OGR in C++: Merge multiple vector datasets to single

texts:
I need to merge multiple kml files into a single one.

I found that I could use "ogr2ogr -append" but I have to do it through c++ and my Google-fu hasn't been successful so far.
So, any ideas?

Also, how to add new fields(and manually values in that field) in a kml file (through C++ again). Probably something to do with OGR_Fld_Create?


texts: Projecting to EPSG 4326, WGS 84 Yielding Different Results in QGIS, ArcMap

texts:
I am projecting a shapefile from "NAD 1983 2011 StatePlane New Hampshire FIPS 2800 Ft US" to "EPSG 4326, WGS 84." However, the result output by QGIS seems to differ by about 3.25 feet from the ArcMap output. The QGIS shapefile is also about 20% smaller in file size. I tried the same function in FME and it matches the QGIS output spatially, but has the same file size as the ArcMap output. I tried in Global Mapper and the output is about .001 feet North of what I get with QGIS and FME

Why don't these match up?

Figure below: Blue is ArcGIS and green is QGIS

Original Vertex: (1082696.828, 1021749.922 feet)

QGIS Vertex: (-71.28406755831366581, 45.30234233581942505)

ArcMAP Vertex:(-71.28406837515765915 45.3023515785411135)


texts: I'm running qgis 1.8, and trying to load a sid file into my map. It imports fine, but seems to ignore the sdw file that comes with it, as the raster shows up without any apparent georeferencing. Is there a way to force the info in the sdw file to be applied to the sid? Thank you for any insight.


texts: ST contains not working properly with trigger

texts:
I have table1, on which people do regular inserts, and table2, that needs to get updated based on what was inserted to table1.

table1 has

    id integer
    type text
    ref text
    geom polygon geometry

table2 has

    idb integer
    status integer
    ref text
    geom point geometry

Users can sometimes do bulk inserts into table1. My trigger function to update table2 is

 CREATE OR REPLACE FUNCTION table2update()
 RETURNS trigger
 LANGUAGE plpgsql
AS $function$

BEGIN
update table2 a set status = 50, ref =  new."ref"
where new.type = 'POLE' and new."ref" is not null and st_contains(new.geom,a.geom);
RETURN new;
END;

$function$
;

trigger

CREATE TRIGGER derp BEFORE INSERT ON table1 FOR EACH ROW EXECUTE PROCEDURE table2update();

People do their inserts on a software called QGIS but it looks like the following

insert into table1(type,ref,geom)
values('POLE', 'derp','geom'),
('not pole', null, 'geom')

Expected behaviour, all the rows in table2 contained by the polygons in table1 will have the ref updated, as well as the status, but only where the type of inserts in table1 are 'POLE'

Actual behaviour, all the rows get updated with the status, but the type is completely ignored - in other words, the ref stays null but the status changes. When a user pastes a single row of type 'POLE', it works as intended however I do not want to encourage people not inputting all the data

How to fix the code?


texts: Print/obtain all vertices in fc1 if list in fc2

texts:
I have a source polyline feature class and a target polyline feature class. Eventually, I want to replace the target feature geometry with the source feature geometry while preserving existing attributes in the target feature. In the code below, I am trying just to print values based on a condition that looks at a fieldname in the target class. The code runs without error but does not output or print anything (I expected it to find both string values in my list). Any ideas?
import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

sourcelayer = arcpy.mapping.Layer('sourcelayer')
targetlayer = arcpy.mapping.Layer('targetlayer')

d1 = {}

with arcpy.da.SearchCursor("sourcelayer", ["OID@","STICK_ID", "SHAPE@"]) as curs:
    #iterate
    stIdCount = 0
    pointCount = 0
    for stId in curs:
        #update dictionary
        #da.[stId[1], [])
        stIdCount += 1
        stickid = str(stId[1])
        print "\n"
        print ("\t"*4), stickid
        for point in stId[2].getPart(0):
            pointCount += 1
            print str(stickid), "    ", point.X, " ", point.Y
            d1.setdefault(key, []).append([point.X, point.Y])
print "\n", "Total: {0} wellsticks and {1} vertices".format(stIdCount, pointCount), "\n"


texts:
I have 250 shapefiles of large area. I have a area shapefile against which I want to clip the small portion of all 250 shapefiles. I tried using transformer Clipper individually on each shapefile as reader and output on individually writer.

As it is lengthy process to assign Clipper to 250 readers individually.

Is there any process to clip all shapes using single clip transformer and output to 250 layers.


texts: Exported layer shows zero geometries using QGIS

texts:
QGIS Desktop application (version X) was used to query a layer and then export result geometries as shapefile. QGIS Desktop application (3.10.12-A Coruña) is used to open the shapefile and view geometries but QGIS shows no features and says that feature count is zero.
Exported shapefile is > 1GB so there must be at least one geometry exported right?


texts:
i try to follow this to Install Third-party Python Modules in QGIS (Windows).my QGIS version is the latest 2.18.10.

steps :

donwload ez_setup.py

2.copy paste at C:\Windows\System32

3.OSGeo4W shell Run as Administrator

4.python ez_setup.py

and i take that message :

any idea why ?

if i try to :


texts: Not all polygons removed when cutting from vector layer in QGIS

texts:
I have created ten polygons representing lakes on a vector layer, and I am trying to cut them from an imported and re-projected vector layer of administrative districts. When I use the 'Cut with polygon from another layer' feature, only nine lakes are removed. The last one overlaps two districts but is not cut from either of them.
Grateful for any help at all.


texts:
I am working on a project to answer the following geographic question: How much have slum areas increased in Mumbai, India, since the 1970's (or earliest data).

I have downloaded LandSat imagery and plan to perform a series of raster classifications to map slum areas. I have imagery from 1975 (LandSat 2), 1993 (LandSat 5), and 2014 (LandSat 8).

Is my approach to this question practical given the data I have access too? I feel I can work with my latest data set, but how much accuracy can I classify images from LandSat 2 and 5? I know that the cell size of LandSat 2 is pretty course(60m). I am using ArcGIS software and would like to know how more experienced GIS professionals would approach this?


texts: Field values missing from quantity symbology options in ArcMap?

texts:
I am using ArcMap.

I joined table data to the latest TIGERLine tract shapefile for my county and exported the data. When I attempt to create a choropleth/thematic map, many of the fields in the table are not visible/options to choose from the Value Field in the Quantities symbology. E.G. I cannot find/choose 2015 Median Income as a value for a graduated color map.


texts:
I created vector polygons from a raster image and now I want to merge adjacent polygons based on a distance. I used the Aggregate Polygons (Cartography) tool after smoothing and running a Repair Geometry tool in ArcGIS.
However, the tool keeps failing. Dissolve tool suggested here doesn't work, because polygons don't have a common attribute. Could someone please help to sort this out or suggest a solution with other software such as QGIS or SAGA.


texts: Why is my DEM layer turning black when I apply multidirectional hillshade effect?

texts:
I'm quite new to QGIS, and am experimenting with the hillshade effects. I've imported a DEM layer, and am simply trying to render it as a multidirectional hillshade. Every time I click the "multidirectional" option and then apply the new settings, the hillshade layer becomes completely black. I've tried changing all of the other settings, and I am able to apply them perfectly fine on their own/until I select the multidirectional option. Any suggestions as to what might be going on? I know there must be a simple answer to this, but I haven't been able to find anything online yet. Thank you in advance for any advice.
Harley


texts: Divide polygon features to reduce feature size

texts:
I am attempting to load FCC Mobile Deployment Form 477 Data into Google BigQuery for analysis and combining with other data I am using. I converted the shapefiles into new-line delimited geojsons using ogr2ogr and concatenated the files for each data release. Unfortunately, BigQuery has a 100MB limit for each row (feature). Initially, I tried to fix this by breaking multipolygons into single polygons, however there are still a large number of features that are much bigger than 100MB (at least one is 500MB). I've attempted simplification at multiple tolerances and the final size is still too large.
What's the best way to divide up the giant polygon features into smaller features that can be loaded into BigQuery?


texts: Determining the ID of the largest area(s), by bounding box(es)

texts:
Here's a curly one!

EDIT: As an alternative the modal raster value inside each bounding box would work, and might possibly be faster. Still unsure how to do it though. Perhaps Zonal statistics using Majority. Then I'll need to filter by the maximum count.

Background Info
Currently using ArcGIS Pro.

I currently have a grid, and a map of vegetation in raster format.
I want to find the most common vegetation by area, in each grid square.
The dataset is probably a couple of hundred thousand rows long - and likely isn't nicely sorted.

Current Process
Vegetation Raster to Vegetation Polygons (Raster to Polygon)
Pairwise Intersect Vegetation Polygons with Grid Polygons to create VegetationGrid_Polygons
Calculate Geometry Attributes on VegetationGrid_Polygons to give Area
Stuck.

Starting Table
Now I should theoretically have a table that looks like:
VegetationID, GridID, Area
1, 10, 2
2, 10, 8 <-- Notice this, there are two VegetationID's #2
2, 10, 4 <--
5, 10, 4
1, 5, 10
3, 5, 46
5, 5, 4

What tool(s) can I use to create a table that has the sum area of each VegetationID by each GridID?
VegetationID, GridID, Area
1, 10, 2
2, 10, 12 <-- This is now 12
5, 10, 4
1, 5, 10
3, 5, 46
5, 5, 4

Then, how to find the VegetationID with the largest area value by GridID?
VegetationID, GridID, Area
2, 10, 12
3, 5, 46


texts: Adding values in two columns to new column on sync table in CARTO?

texts:
We have field data coming in from an app platform and we want to drive a CartoDB visualisation to help with planning.

A bunch of rangers go out surveying sites to count crocodiles. For each site, the app collects:
- Number of Male crocodiles (num_croc_male)
- Number of Female crocodiles (num_croc_female)

We don't actually want our people standing there in the dark doing sums, so they upload the data in this format straight to the cloud, and this is then POSTed over to our CartoDB.

In the CartoDB visualisation, we want to show the two cols above, and we also want: Total Number of Crocodiles (num_croc_sitetotal).

Since this integer will also be used for styling the display of each site, I want the value stored in a column.

And since the data is being refreshed after each patrol, I need this to be a stable transform that persists through new rows being added.

It seems like I can't find the tutorial or walk through for how this might be done...and my SQL is nearly non-existent. I was sort of hoping CartoDB had some sort of formula-building interface that would shield me from the full blast of SQL syntax.

Could anyone provide an outline of the solution steps I should be following here?


texts: Saving selection as new layer in QGIS 2.4?

texts:
I've selected a number of polygons from one layer and want to create a new layer. Because I can't right click using a mac how do I find "Save Selection As". Its not in the list in the "Layer" drop-down from the top toolbar.

I'm using QGIS 2.4 with Mac OS X 10.9.4.


texts: How to speed up a complex process of compositing multispectral images in python?

texts:
I would like to speed up my process.
I have a one month time series which consists of 5 images (Sentinel-2) taken over several days. Each image is made up of the same bands that I selected. After several steps I've transformed each image into into a dataframe and, after, I've concatenated all in a multindex DataFrame. Here is the result.

(10 bands x 6 dates)
(dates = %y-%m-%d)
Each column of the first level is a date, each column of the second level is a band, each row is a pixel.
Now it gets complicated.....
I built a code that, for each pixel, I find the three dates in which the values of the B8A band are lowest (practically the three lowest values). Next, between these three, I choose the one with the highest B12 band. Thus, the date in which the B12 is highest will call all the band values of that specific pixel in the new composite DataFrame (since is a multindex, selecting the date brings with it not only the value of the B12 but also the values of the other bands connected to that date).
In the meantime, I list all the dates for each pixel in exact order, so that I can associate them for future analysis.
Here the code:
columns = ['B2', 'B3', 'B4','B5','B6','B7','B8','B8A','B11','B12']
im_04Jul = pd.DataFrame(image1, columns=columns)
im_09Jul = pd.DataFrame(image2, columns=columns)
im_14Jul = pd.DataFrame(image3, columns=columns)
im_19Jul = pd.DataFrame(image4, columns=columns)
im_24Jul = pd.DataFrame(image5, columns=columns)
im_29Jul = pd.DataFrame(image6, columns=columns)

# dates = DateTimeSerie  (first level)
df_Jul = pd.concat([im_04Jul, im_09Jul, im_14Jul, im_19Jul, im_24Jul, im_29Jul], axis=1, keys=dates)

#creating the DataFrame with only one band
def conc_band(band):
    return pd.concat([im_04Jul[band], im_09Jul[band], im_14Jul[band], im_19Jul[band], im_24Jul[band], im_29Jul[band]], axis=1, keys=dates)

B8A = conc_band('B8A')
B12 = conc_band('B12')

#definition of the max B12 rule
#return the dates (first level) where the B12 is higher
def max_swir(B12, cont):
    a =B12[cont:cont+1][date_ultima].item() #.item() take only the number
    b=B12[cont:cont+1][date_penultima].item()
    c=B12[cont:cont+1][date_terzultima].item()
    if a > b & a > c:  #For pandas "and" and "or" are considered ambiguous so you should use "bitwise" "|" (or) or & (and) operations:
        return date_ultima
    elif b>a & b>c:
        return date_penultima
    else:
        return date_terzultima

#start the composite construction
composite = pd.DataFrame()
dates_composite_list = []
cont = 0
while cont < len(df_before):
    df = B8A[cont:cont+1].sort_values(axis=1, by=cont, ascending=False) #sort the minimum value at the last
    df = df[df.columns[-3:]] #select the last 3 values
    #obtaining the name of the three last dates to use them to select the other bands
    date_ultima = df.columns[-1] #the name of the first last
    date_penultima = df.columns[-2] #the name of the second last
    date_terzultima = df.columns[-3] #the name of the third last
    date = max_swir(B12, cont)
    new= df_Jul[cont:cont+1][date]
    composite=composite.append(new)#HERE THE COMPOSITE!!!!!
    dates_composite_list.append(date) #here we create a list with the respective date for each pixel
    cont+=1

The code, although theoretically complex, works very well. The problem is that it is very slow since the pixels are more than 30 million per image and then I have to apply it over several months. Is there a faster way?


texts: Having trouble publishing from map documents

texts:
I am attempting to create an operational layer for the Flex viewer. I am having difficulty getting my mxd to publish correctly (not visible when on server). Here is the error message I get after I create a layer package and try to publish using publish tool in ArcMap 10.

  Possible cause of connection failure: the soc cannot resolve the path
  to the data, or the soc account has not been granted read permissions
  of the data. In the case of ArcSDE, direct connect may not be
  configured properly, or the soc account may not have permission to
  access the database.

Any one know what I should do?


texts: Using forceToLineString() with GDAL and Python?

texts:
For my script I'm looking for a method to convert a MultiLineString into a simple LineString.

Example:

MULTILINESTRING ((0 0,0 1),(0 1,0 2),(0 2,0 3,0 4))

==>

LINESTRING (0 0,0 1,0 2,0 3,0 4)

I found a GDAL-method called "forceToLineString" but I don't know how to use it.
The way I tried does not work:

from osgeo import ogr
multiline = ogr.Geometry(ogr.wkbMultiLineString)

line1 = ogr.Geometry(ogr.wkbLineString)
line1.AddPoint_2D(0,0)
line1.AddPoint_2D(0,1)

line2 = ogr.Geometry(ogr.wkbLineString)
line2.AddPoint_2D(0,1)
line2.AddPoint_2D(0,2)

line3 = ogr.Geometry(ogr.wkbLineString)
line3.AddPoint_2D(0,2)
line3.AddPoint_2D(0,3)
line3.AddPoint_2D(0,4)

multiline.AddGeometry(line1)
multiline.AddGeometry(line2)
multiline.AddGeometry(line3)

multiline.forceToLineString()

print multiline

Can someone help me?


texts: Calculate nearest feature with a barrier

texts:
I have three shapefiles: a) polyline with rivers and streams (also as polygon available to me), b) polylines of highways and roads and c) point file with locations in the riparian vegetation. I want to calculate the distance of each point to the nearest road or highway,
on the same side of the river.

This means, that the tool "NEAR" in ArcGIS is not appropriate as it calculates the distance to the nearest road regardless of the side. What I am looking for is to incorporate a barrier, or anything similar in its outcome, that only the distances from the locations to the roads on the same side of the rivers are calculated.

What would you suggest?

I am using ArcMap 10.0.


texts: Convincing organization to move towards Open Source GIS?

texts:
After spending the better part of the year on this site I am seeing how big the open source GIS community is. I used to just think it was ESRI and nothing else. That said, how do you begin to nudge or convince others that open source programs might be more beneficial than ESRI products?

I'm not specifically asking about features in one vs the other, but actually trying to convince higher ups that a move to a similar product that would provide comparable results in terms of productivity, efficiency and end product.

I guess what it boils down to in the end, is, at the end of the day, does the cost savings outweigh the time that it might take to get GIS staff acclimated to open source software?


texts: Creating buffer around stream network within raster layer?

texts:
I have a raster layer of a stream network with two values (1 no stream, 0 stream). I need to create a 100m buffer around the entire stream network that lies within a specific watershed.shp. I have access to Spatial Analyst.
I've tried the expand tool and the polygon to raster tool with zero success.


texts: Extracting city information only from OpenStreetMap

texts:
I need to extract City as well as hierarchy of admin entities for selected countries.

My understanding for now is that I'll have to get node  primarily and find a clever filtering based on tags. And possibly include some ways as they are often used to represent administrative boundaries.

I also need to store these info in a specific database design in Postgresql.

Would osm2pgsql be the most adapted tool to filter out information and directly import in my db design?

Or would overpass be necessary anyway to do the filtering because I feel it'd would be a much more powerful tool in this regards?

I'd then need to program import to DB based on received data, but I'm not even sure osm2pgsql would be flexible enough to import in my design.


texts: Exporting raster as .tif with "use renderer" and "force RGB" produces inaccurate colors?

texts:
In ArcMap 10.2 I have a DEM and my aim is to change its color ramp using ArcMap, then export this colored raster as a geotiff (8 bit).

I checked "use renderer" and "force RGB" options but the results are not good (EDIT: when viewed in ArcMap). With some color ramps, the results are really off, and with others, they are slightly off.

For example:

Is this a bug, or is this an artifact stemming from the the original image being simplified into an 8 bit image? (Does ArcMap stores its color ramps with more than 8 bits?)

it seems that results only look strange in ArcGIS - if viewed in other applications, colors look right.

So why is ArcMap displaying the exported rasters inaccuractly?

Do I need to change some settings somewhere to have them display correctly?


texts: How to configure mapproxy and mapnik for getfeatureinfo requests

texts:
I'm wondering if its possible to configure mapproxy and mapnik for getfeatureinfo requests. I've made this combination work for getmap requests, but when I request getfeatureinfo I receive layer xxx is not queryable. I've tried setting featureinfo:true, but no luck there.

Any help appreciated.


texts: I have many "old" (10+ year) MapInfo workspaces that I want to reopen but cannot locate some of the tables. MapInfo will not open the workspaces without these tables and the data/objects I want only exist in the editable layer of these old workspaces. Is there a way of overiding MapInfo's inability to open the workspace without these tables??


texts: Capturing specific error instead of arcpy.ExecutionError?

texts:
How can I specify this script to capture specifically the error 001566, instead of the current general execution error?

  Error 001566:Failed to create raster dataset. The reason could be: The
  clip feature is outside the raster extent

try:
    arcpy.Clip_management(lyr,rectangle='#',out_raster=out_feature, in_template_dataset = clip_feature, nodata_value="256", clipping_geometry="ClippingGeometry", maintain_clipping_extent="NO_MAINTAIN_EXTENT")
except arcpy.ExecutionError: