texts: ArcGIS: fishnet with rotation angle and coordinates error

texts:
I have a micro plots map shapefile with NDVI data. I want to create a fishnet of polygons to divide the map into the micro plots. I have being able to create a fishnet with polygons but I cannot place it on top of my map.

I know my map has a -32 angle rotation. I have tried to modify the coordinates but still get the wrong placement for fishnet as you can see in the image.

I am now using ArcMap 10.5


texts: Convert height Geoid EGM96 to ellipsoidal?

texts:
I'd like to convert the height data which is in EGM96 to Orthometric projection.

Perhaps there are tools, like ArgGIS (I cannot find if ArcGIS can do this), or any formula that may help me.


texts:
Is there a shapefile anywhere showing Mexican boundaries (national & state) that are NOT clipped to the coastline? Either polygons showing the maritime boundary, or polylines showing only the land borders without the coastline would work.

I prefer to use the INEGI vectors, but I've only found those with coastlines.


texts: I was wondering if there is a website that showcases Openlayers and PostGIS 3D in Action


texts: Changing the default GeoNetwork Map Viewer basemap (OSM) Center

texts:
By default, GeoNetwork Map Viewer renders the base OSM map at center 0 with zoom 2. I would like to change the center to a different point (like Paris) with a zoom of 6.

How to accomplish it by changing the 3.8.1 geonetwork core source directly?

From mapsManager.js, I tried changing center and zoom to different value under:

createMap: function(type) {
          var config = gnMap.getMapConfig()['map-' + type];
          var map = new ol.Map({
            layers: [],
            view: new ol.View({
              center: [0, 0],
              projection: gnMap.getMapConfig().projection,
              zoom: 2
            }),

I rebuilt the war and deployed it to app server, but it does not seem working.

Also saw the post on:
Focusing GeoNetwork default map on Canada?. But not clear what exactly change to the config-viewer.xml

Any insights?


texts: Initialising QgsFeature with predefined field names using PyQGIS?

texts:
I need to create function that assembles QgsFeature from scratch: takes id, dictionary of attributes, and geometry as input. The issue is that feature that is being created doesn't belong to any layer and not going to before heavy processing that needs this feature to have attributes callable by name. Attributes are fixed: 2 fields, one is string and other is real. From documentation I can see that QgsFields() can help me to assign fields to the feature being created. Though I don't get how does QgsFields() work.

#att_dic = {u'value': value, u'zone': str(some_data)}
def makeFeature(fid, att_dic, geom):
    feat = QgsFeature(**some QgsFields(.....) magic??? **, id=fid)
    feat.setGeometry(geom)
    for name, value in att_dic.iteritems():
        feat.setAttribute(name, value)
    return feat


texts: Loading layer styles; (*.qml file) is there a shortcut or toolbar icon

texts:
Currently, to load a saved style you have to open the layer properties, scroll to the bottom and click on the style tab. As this is a regular thing I need to do when creating new maps I was wondering if I am missing and easier or quicker way to do this?
I have the layer styling window open all the time, but this window doesn't have access to the style button I need to load specific styles I regularly need when starting a new map.


texts: Enforce that all "Enforce expression constraints" of child tables are valid before I can save the parent feature

texts:
In QGIS 3.16 I set up a relation between two tables.
In the child-table I set a constraint-expression and checked "Enforce expression constraints".
If I add one or more features to the child-table which have attributes, which are not valid, QGIS recognizes them as "invalid" as expected but:

I can still toggle editing mode for child layer - although some data in child-layer is "invalid"
I can still save the parent-layer ("OK"-Button is not disabled)  - although some data in child-layer is "invalid"

How can I prevent the user from saving invalid data in both cases without python-scripts (like How can i prevent a form from closing unless data is accepted by validation) ?
(I thought to set up a constraint-expression and "Enforce expression constraints" in the parent-table, which checks all child-feature data - but there seem not to be a way to do this.)
I'd like not to use Python because:

I have a parent-/child-relation (which is -as far as I know - very complicated to realize with custom forms) and
I'd like to use this project in QFIELD, which has - as far as I know - limitations with custom-forms - so I use drag and drop-forms


texts: ArcPad 10 Studio, from form to vbscript back to form

texts:
There is very little documentation on using vbscript within a form. Any help with syntax for ArcPad 10 would be appreciated.

Within my subroutine, I want to grab some numerical data entered in a control box on the form. Will this Set command grab the data that was input and assign it to the variable el1?

Sub getdata
Dim el1
Set el1 = EDITFORM.Pages.Item("page4").Controls.Item("element1")
End Sub

To bring data from the vbscript back to a control box on the form, assuming I have calculations done correctly for "scoreA"

    EDITFORM.Pages.Item("Results").Controls.Item("sumEl").Additem scoreA, scoreA

Getting very close.  Just need to get the data from the script back to the form. The above command line gives a runtime error "Class Doesn't Support Automation"


texts: Convert Lat/Long to UTM specifying which zone to use (C++)

texts:
I need to convert Latitude/Longitude coordinates into UTM, in C++, specifying which zone to use.

I am working with locations that fall in two zones, and would like to project them all into one zone - since they are not too far from the border, I am hoping the error introduced by doing so won't be too much.

So far I had been using Chuck Gantz's library, which converts LL to UTM, but it does not support conversion from LL to UTM specifying which zone to use.


texts: It possible to query/convert LONG RAW (SDEBINARY) from SQL?

texts:
Curently we have ArcSDE 9.3/Oracle10, and layers in SDEBINARY.

It is possible to make queries with SDE.ST_GEOMETRY or SDO_GEOMETRY commands ?


texts: Add Geometry Attributes doesn't work for feature classes but does work for shapefiles

texts:
I'm trying to add start and end point coordinates to a line feature using the LINE_START_MID_END geometry property of the Add Geometry Attributes tool. Every time I use the tool on a feature class in a geodatabase I get an error message:

<Error>

Traceback (most recent call last):

  File &quot;c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\AddGeometryAttributes.py&quot;, line 292, in &lt;module&gt;
    AddGeometryAttributes(fc, geomProperties, lUnit, aUnit, cs)
  File &quot;c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\AddGeometryAttributes.py&quot;, line 49, in AddGeometryAttributes
    fields = CreateOutputFields(fc, geomProperties, hasZ, hasM)
  File &quot;c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\AddGeometryAttributes.py&quot;, line 244, in CreateOutputFields
    arcpy.da.ExtendTable(fc, &quot;OID@&quot;, narray, &quot;_ID&quot;)
RuntimeError: General function failure [FeatureClassInQuestion]
</Error>

<Error>Failed to execute (AddGeometryAttributes).</Error>

</Errors>

when I convert the feature class to a shapefile and run the Add Geometry Attributes tool on it, it runs fine and produces the points.

I've tried disabling M and Z values in environments without success. I suppose I could just convert the feature class to a shapefile and convert it back to solve this issue but I'm curious to see if anybody else has run into this and found a simpler way to fix it?


texts: Published feature service with massive features drawing too slowly

texts:
I have published a feature class with a massive number of polygons with large vertices on ArcGIS Online. When opened and viewed in Map Viewer, it drew acceptably slowly.

I noticed it was hosted published as a Feature Service. I opened the Developer tool on the web browser. Every zoom in/out and pan triggered a massive data (10MB plus) transfered from the FeatuerServer on ArcGIS Online. What was returned was a large collection of feature geometry. That's why it took ages to draw on the map.

I am wondering if the hosted feature class can be published as Map Service or WMS so that the data transfered could be simply images like what is returned by a "export" capability for a MapServer.

By the way this hosted feature layer woule be updated quite often, e.g. twice a day.


texts: Accessing all layers displayed on map in Mapbox

texts:
I want to access names of all layers that I display on the map using Mapbox GL with React.
Is there any specific method that enables me to access that?


texts: Road Names - Multiple Lines and Single Lines

texts:
Fairly new to the GIS world. I have been creating maps for my employer using OS Vector map data. Up until recently everything has been going smoothly - I only use vector based data so that I can adjust the colours and styles of the map. However, the only issue I have is the road names when compared to a rasterized version.

The rasterized version (Above) includes all road names using a variety of different placements. How can I replicate this using vectors and QGIS?


texts: How to get decimal degrees instead meters from X Y coordinates in QGIS?

texts:
I created a new project with decimal degrees system. Then I created a new layer and created a point on it. When I click on my point I see in attributes X and Y attributes in decimal degrees but when I try to create a new field with function concat($y,' ', $x ) it displays value in meters instead degrees. What I have to do for get value in decimal degrees?
QGIS 3.12


texts: Seeking open raw LiDAR or laser scanner data?

texts:
I'm looking for aerial LiDAR or terrestrial laser scanner data in laz or las format preferably from an FTP site.

Does anyone know which countries have that data available to the public?

From my research it doesn't seem to be easy to come by.


texts: Plotting square grids on a map and extracting each grid information using R

texts:
I have a data set which looks like this

LATITUDE LONGITUDE VALUES
28.26734 76.82973 -56.7879
28.26644 76.82972 -58.22462
28.26508 76.82969 -58.4211
28.26778 76.83076 -55.75333
28.26733 76.83075 -58.55153
28.26507 76.83071 -56.38619
28.26912 76.83129 -56.11011
28.26732 76.83126 -58.17415
28.26687 76.83125 -59.77212

I want to plot the "Values" on the map wrt its "LAT" and "LONG" and create square grids on that map. When the grid is created, each grid will contain some rows of the data set i.e. the "Values" which are plotted. and I then want to extract the values each grid is having. So the fourth column in the dataset would be "GRID" and each row belonging to a particular grid...Of course multiple rows can take single grid values.

I am able to plot these points on the map but not able to create the grid.


texts: Changing the Location of Python Expressions QGIS3

texts:
I've recently started using QGIS and I've written a few custom functions to use with expressions in the Function Editor (see image below)

From what I can find these appear to be saved in:

C:\Users\username\AppData\Roaming\QGIS\QGIS3\profiles\default\python\expressions

I would like to share these with my colleagues. Is it possible to simply change the location these custom expression functions are read from? Or add another location they could be read from?

I would prefer the latter so that everyone could have access to the same general custom functions while still being able to access/write their own custom functions without everyone else having access to them.


texts: Adding a column and update values using SQL in QGIS

texts:
I have a vector file "data.shp" containing the following fields:

ID, value1, AREA

I want to add a new column "value2" and set it as value1/AREA. Since I want to do that for many columns and files I want to use SQL.

I tried with virtual layer in QGIS with no success. I wrote in the query the following :

ALTER TABLE data
ADD value2 float;
UPDATE data
SET value2 = value1/AREA

Are alter and update commands working with virtual layer? Is there any other way?


texts: ArcGIS Toolbox, parameters for python script not as positional arguments

texts:
I have an ArcGIS Toolbox with a bunch of parameters. A python script is used to perform some raster operations, based on this parameters. I am using Python argparse in the script to get the parameters. I do this because the program was originally developed as a command line tool using ArcPy. With argparse I can show some useful help text and use command line switches where the position of the argument in the argument list does not matter. Of course argparse offers much more advantages.

ArcGIS is passing the parameters as positional arguments. So I can't use parameters like -k or --doSomeThing=100

Using a checkbox for example, that is ticked, this will pass a "true" to my script but I would like to have a "real" command line switch like "-k" or something like that. Has anyone an idea if this is possible? Maybe using tool validation?


texts: Trigger to SELECT schema rows based on WHERE condition

texts:
I am learning about triggers.

Generally it seems triggers are employed on the backend to perform some action on a schema. I was interested if it was possible for a trigger to return data on a SELECT statement to the main application.

In the table I'm working with we have several string columns followed with a geom (geometry::polygon) column which would contain the WKBT LAT/LONG values.

I'd like to return the rows which meet some WHERE clause AND also return the geom [LAT/LON] in WKBT format. By return I am interested in returning data to the main application.

I seek any examples.

Currently I build a list of all our column names {colA, colB, colC, colD} and then append the following as the last entry: public.ST_ASText(geom).

So I obtain a SELECT as follows....
SELECT colA, colB, colC colD, public.ST_AsText(geom) FROM schemaA WHERE .

So the crux of the problem here I cannot use SELECT * since I must convert the geom (LAT/LON pairs) to WKBT and this must be in the select along with every other column name defined, whether I use that column data or not.  I definitely require the geom column in WKBT format.

This is passed via our DEVART interface.
I was wondering if there is a better way to do this,


texts: QGIS - Zero values after joining .shp to .csv

texts:
I'm using Australian Bureau of Statistics (ABS) data to create thematic maps. Everything works fine when using POA (post code) data. I've created a .csvt sidecar file for the CSV. After the join I can make graduated and categorised maps fine.

I follow the same procedure for SSC data (suburbs) but after the join QGIS returns values of 0 for all fields when trying to make graduated maps. I use QGIS 3.x

Posctcode Screenshots:
Fields:Postcode fields
Attribute Table:

Produces Map:

But for SSC:
Fields:

Attribute Table:

Produces Zero Values:


texts: Attribute Table Subtract Value From Previous Value

texts:
I have a river shapefile and need to ensure every vertice has a lower elevation that the previous uphill vertice.

I have converted the polylineZ to points using ET GeoWizards.

This gives me a Z value for every vertice.

I have added a field called ElevDiff.

I need to calculate the difference in elevation values from the previous vertice to the next in the attribute table.

If the value is negative, I know the line does not flow downhill at this point.

This is an easy task in Excel but I am wondering if there is a python script or easy way to calculate this using ArcGIS.

I do not have any extensions (i.e. Spatial Analyst, 3D Analyst, Newtwork Analyst).


texts: List of all features inside a buffer (Web map)

texts:
I have web map, created in the combination of: ExtJS 5, Openlayers 2, Geoext 2, and Geoserver.

My layers are in postgis linked to geoserver.

My scenario are 1 base map, an 10 normal layers. Actually i have a tree panel and i can view correct all the layers in the viewport.

The idea is the user click a point in the map and in someplace also input the meters to make a buffer of that point. Next i need to get a list of all the features of my 10 layers that are inside that buffer (ring).

How can I achieve this?


texts: Error Installing sample datasets with QGIS 3.4

texts:
Recently tried installing QGIS 3.4 on Windows 10 and somehow it always throws error while trying to install the sample training data sets.

While QGIS itself installed without issues, the data sets failed to load.

Any suggestions?


texts: Printing large number of Points in ArcGIS

texts:
I have bathymetric data which have thousands of points ( exactly 71 K).

Now when I try to print it in a large map size in ARCMAP 10.0 , the layout view freezes and the system is hung , possibly due to a memory issue.

What is the best way to do this, rather than generalizing the points ? my preferred paper size is A0.

Is there a better software to do this ?


texts:
When I use the following code to perform some raster calculations with gdal_calc.py:

gdal_calc.py -A merged.tif --outfile= result.tif --cal="(A)/10000"

I get the following errors:

Module compiled against API version 9 but this version of numpy is 6.
ImportError : numpy.core.multierror failed to import

Python Version :Python 2.7.2


texts: Moving data from columns in a table to rows by a key field?

texts:
I would like to shift data found in a column to a row by Key field. I have a table, within a geodatabase, set up similar to what's below.

KEY_    Value
123     Yellow
123     Blue
456     Red
456     Yellow
789     Green
246     Green
246     Orange

I want to create a table that moves the data in the Value field to rows based on the KEY_ field. The result will look like what's below.

KEY_    Value001    Value002
123     Yellow      Blue
456     Red         Yellow
789     Green
246     Green       Orange

Does anyone have an idea of how this can be done? I tried using the Transpose tool Esri offers but it didn't produce the result I'm looking for. I have basic understanding of Python.


texts: Click on external div link to Zoom to feature (or lat/long)

texts:
I want to click on a link in a navigation div in order to zoom to a location within my leaflet map.

Html:

First link

Second link

Javascript:

//initialize map
var map = L.map('map', {
  zoomControl: true
}).setView([40.7451, -73.9459], 16);

var osm = new L.TileLayer('http://{s}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png', {
  maxZoom: 19
}).addTo(map);

//click and zoom
$('story').on('click', function(){
    // parse lat and lng from the divs data attribute
    var latlng = $(this).data().point.split(',');
    var lat = latlng[0];
    var lng = latlng[1];
    var zoom = 10;

    // set the view
    map.setView([lat, lng], zoom);
  })

I am using:
jQuery v3.2.1
leaflet v0.7.7

When I click on the div nothing happens. Forgive me if there is an obvious solution, I have brand-new to java-script.

I'm attempting to use a solution originally posted here:
Leaflet.js map - Pan too and Zoom in on point


texts:
Is there a way to save the current shapefile, raster, etc selections in QGIS and then bring them up all at once later? For example, I find myself re-creating maps over and over because I need different shapefiles or labels or colors. I would like to save the current extent and selection I have so I can quickly reference back to a map I've already made.

Is my question clear? I have tried asking The Google but can't seem to phrase it right.


texts: How to reset the id number of newly generated PostGIS topology?

texts:
I'm trying to create a clean, new topology for testing by discarding the old topology first, i.e.:

SELECT DropTopology('ch13_topology');

SELECT CreateTopology('ch13_topology', 32610, 0.05);

However, one issue I have is that topology.CreateTopology is giving a different (increasing) id number each time for the topology. The first time, it returns 1, the second time it returns 2, and so on, even though the old topology has been dropped.

I understand it's safe to always increase the id number. But,

Is there a way to make topology.CreateTopology return the smallest id number possible?

In the case shown above, this would be always 1, or whatever number it is at the first call of {drop,create}topology(). I need this so that I don't have to change the id numbers in subsequent SQL statements.


texts: How to query objects within 1000m epsg4326?

texts:
I have a PostGIS table having a column "geom" having a type geometry(MultiPolygon) and the projection epsg:4326 (degree decimal), I want to query all objects within 1000m and I use this query:

SELECT ST_AsGeoJSON(geom) FROM table
WHERE ST_DWithin(geom, 'POINT(long lat)', 1000.0);

but I always ended up querying all objects.


texts: Does power of GPU affect performance of analyses by ArcGIS Desktop?

texts:
Does the power of GPU affect the performance of analyses such as Segment Mean Shift, conversion to/from raster, etc?

Is the GPU used only for caching GUI elements to give a smoother user interface or is it used for processing as well?


texts:
I'm running Geoserver over HTTPS and when I preview a layer with openlayers3 format the page throws an error.

The error point that some resources

i.e. http://mydomain.com/geoserver/openlayers3/ol.js

are blocked because are requested over HTTP.

I tried to change it manually but I can't find the file.

Is there some configuration file I need to change so that the page request the content over HTTPS?


texts: Legend of a raster map with categorical data

texts:
I would like to plot a raster containing 4 different values (1) with a categorical text legend describing the categories such as 2 but with colour boxes.

I've tried using legend such as :

legend( 1,-20,legend = c("land","ocean/lake", "rivers","water bodies"))

but I don't know how to associate one value to the displayed color. Is there a way to retrieve the colour displayed with 'plot' and to use it in the legend?

Initial raster

Raster with legend


texts: Display a layer from geoserver using OpenLayers 3

texts:
I'm wondering how to display a Layer in OpenLayers 3 which I host in Geoserver. I imported the Layer from a PostGIS-DB and it is possible to display it using the geoserver preview option. The Layer uses EPSG:3857.

    <link rel="stylesheet" href="http://openlayers.org/en/v3.1.1/css/ol.css" type="text/css">

    <title>test</title>

The Layer should be displayed in this area: 52.025151, 13.069072


texts: Passing Selection from Combo Box into Definition Query using Python AddIn?

texts:
I am trying to pass the users selection into a definition query in an addin combo box python script.

Basically, I have a feature class with a field called "SITE". I build a python list from that field (using the search cursor and appending the results) and pass it into the combo box items list. The user will choose a value in the combo box box from that list. I want to further pass the users selection into a definition query.

You can see in my script below, I just tried to pass the selection into the definition query, but it crashes the whole combo box and toolbar. I've tried converting it to string, placing quotes around it, etc and no luck. Obviously I am doing something wrong with. I believe it's how I'm passing the selection. As suggestions?

import arcpy
import pythonaddins

sitename = []
for row in arcpy.SearchCursor("Z:\Geodatabases\Test.gdb\Tesfc"):
sitename.append(row.SITE)

class ComboBoxClass(object):
    """Implementation for SolsticeFigureUpdater_addin.combobox (ComboBox)"""
    def __init__(self):
        self.items = sorted(sitename)
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWWWWWWWWWW'
        self.width = 'WWWWWWWW'
    def onSelChange(self, selection):
        mxd = arcpy.mapping.MapDocument("CURRENT")

        for df in arcpy.mapping.ListDataFrames(mxd):
            if df.name == "DFMAIN":
                main = arcpy.mapping.ListDataFrames(mxd, "DFMAIN")[0]
            if df.name == "DFINSET":
                dfinset = arcpy.mapping.ListDataFrames(mxd, "DFINSET")[0]

        for lyr in arcpy.mapping.ListLayers(mxd, "", dfmain):
            if lyr.name == "PDA Boundary":
                lyr.definitionQuery = '"SITE" =' selection

        arcpy.RefreshActiveView()
        arcpy.RefreshTOC()

    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass


texts:
From what I've understood ol.Object classes have observable properties which are detailed here: http://openlayers.org/en/v3.8.1/apidoc/ol.Object.html

This makes sense to me, however what is not detailed here is how to dispatch custom events on classes that extend the ol.Object.

In my current application I create controls that extend the ol.Control class and by definition are ol.Object's. So I can listen to changes to properties that I have set using ol.Object's get/set methods. I've also been able to see how events are dispatched by doing the following for example:

objectOne.dispatchEvent('fired');

I could then listen to his event and run a function accordingly:

objectOne.on('fired', function(e) { alert('fired!'); });

In my code at present I am performing these operations and everything is ok, providing I use the debug version of OL3. As soon as I change to the compiled version dispatchEvent is no longer available.

Why is this method hidden in the production version? Should I be using a different method to fire my own custom events? It seems strange at present that you can listen for events and property changes but you cannot fire your own custom events when using the compiled version of OL3


texts: QGIS GRASS maximum distance to a given feature

texts:
I want to calculate the maximum distance between a point and the boundary of a polygon (representing an animal's home range).

Here is an example of the point layer, always within the polygon layer:

When I needed the minimum distance, I used "QGIS GRASS > v.distance > minimum distance to nearest feature".

Does anyone know how to calculte the maximum distance in QGIS/GRASS?


texts: I have digitised glaciers in Southern Chile and at one stage, my shapefiles were lining up perfectly and somewhere along the way I have obviously changed something and now the shapefiles do not line up with the image. I have tried changing the projections and even when they are set to the same projected coordinate system and they end up nowhere near each other. I can set the image to WGS_1984_UTM_Zone_19N and the shapefile to WGS_1984_UTM_Zone_24N, the polygons are on the image but they are just slightly off to the right of where they should be.


texts: How to add a layer file to an ArcGIS map document programmatically?

texts:
I want to create a feature layer from a layer (.lyr) file and add it to my ArcGIS map document programmatically. I am using C# and ArcGIS.

I have a layer file and a shapefile.


texts: addLayer with switchLayer in OpenLayers 3?

texts:
I found this example for switching layers in OpenLayers 3.

function switchLayer()
 {
  var checkedLayer = $('#layerswitcher input[name=layer]:checked').val();
  for (i = 0, ii = layers.length; i < ii; ++i) layers[i].setVisible(i==checkedLayer);
 }

$(function() { switchLayer() } );
$("#layerswitcher input[name=layer]").change(function() { switchLayer() } );

var layers = [];
layers[0] = new ol.layer.Tile({ source: new ol.source.MapQuest({layer: 'sat'}) });
layers[1] = new ol.layer.Group({ layers: [ new ol.layer.Tile({ source: new ol.source.MapQuest({layer: 'sat'}) }), new ol.layer.Tile({ source: new ol.source.MapQuest({layer: 'hyb'}) }) ] });
layers[2] = new ol.layer.Tile({ source: new ol.source.MapQuest({layer: 'osm'}) });
layers[3] = new ol.layer.Tile({ source: new ol.source.OSM() });

var map = new ol.Map({
  target: 'map',
  controls: ol.control.defaults().extend([ new ol.control.ScaleLine({ units:'metric' }) ]),
  layers: layers,
  view: new ol.View({
    center: ol.proj.transform([0, 10], 'EPSG:4326', 'EPSG:3857'),
    zoom: 3
  })
});
</script>

<input type="radio" name="layer" value="0" checked> MapQuest Satellite

<input type="radio" name="layer" value="1"> MapQuest Hybrid

<input type="radio" name="layer" value="2"> MapQuest OSM

<input type="radio" name="layer" value="3"> OSM

This is great for having multiple base layer maps...  In addition to this, I also need the ability to add layers (addLayer) to this same map that allow me to add multiple feature/vector layers to these maps via checkbox section.  Ideas?
I was thinking something like this:

$("#layerswitcher input[name=layer]").change(function() { addLayer() } );

var layers = [];
layers[4] = new ol.layer.Image({
source: new ol.source.ImageWMS({
    url: 'http://demo.opengeo.org/geoserver/wms',
  params: {'LAYERS': 'topp:states'},
  serverType: 'geoserver'
      })
    })
  ];
layers[5] = new ol.layer.Image({
  source: new ol.source.ImageWMS({
    url: 'http://demo.opengeo.org/geoserver/wms',
    params: {'LAYERS': 'countries'},
    serverType: 'geoserver'
    })
  })
];
layers[6] =   new ol.layer.Image({
  source: new ol.source.ImageWMS({
    url: 'http://demo.opengeo.org/geoserver/wms',
    params: {'LAYERS': 'point_clouds'},
    serverType: 'geoserver'
   })
  })
];

but it isn't working.


texts: Aligning many rasters using PyQGIS or Python?

texts:
I want to align many rasters in some pixel size, extent and projection system using Python or PyQGIS.

First think is to use GDAL:

gdalinfo (to find information from reference raster)
gdalwarp (to convert projection system and pixel size)
gdal_translate (to convert extent)

They work, but it is not easy to use those tools together in Python and they need much time to finish this work. Also, it can be done easily using QGIS and Align Rasters Tool.

Does a tool doing this work using PyQGIS or Python exist? (I want to work more programming automatically)

UPDATE

I find this code from this question:

from osgeo import gdal, gdalconst

inputfile = #Path to input file
input = gdal.Open(inputfile, gdalconst.GA_ReadOnly)
inputProj = input.GetProjection()
inputTrans = input.GetGeoTransform()

referencefile = #Path to reference file
reference = gdal.Open(referencefile, gdalconst.GAReadOnly)
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)
x = reference.RasterXSize
y = reference.RasterYSize

outputfile = #Path to output file
driver= gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile, x, y, 1, bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)

gdal.ReprojectImage(input, output, inputProj, referenceProj, gdalconst.GRA_Bilinear)

del output

That code works fine except raster clip. Any idea how to update this code to clip input raster to extent of reference raster in code?


texts:
I want to elevate a certain area of terrain by 126 meters. I've created a polygonal shapefile layer with a decimal field titled "elevation" set to 126. My plan was to rasterize the polygon using GDAL Rasterize, set the null areas to 0 using GRASS r.null, and then add the result to the original elevation raster. When I used GDAL Rasterize, the resulting raster was transparent, I assume because it was constituted completely by null values.

How can I turn my polygon into something that I can add to the original elevation raster?


texts: Escaping spaces in paths for use in batch file

texts:
I am making a MapInfo program in MapBasic. I am calling OGR2OGR in order to convert some files. To do this I am generating a batch file and then running it. My code looks something like this:

RunStr = "OGR2OGR -f ""ESRI Shapefile"" """ & ExportFolder & "SHP\"" """ & TableFileName & """ -overwrite -skiperrors"
SHPMaker = ExportFolder & "SHPMaker.bat"

Open File SHPMaker
  For Output
  Access Write
  As #1
Print #1, RunStr
Close File #1

Run Program SHPMaker

My problem is that ExportFolder and TableFileName can potentially include spaces. My questions is how should I go about escaping spaces? I have tried both adding carets (^) before the spaces and putting double quotes around the paths -neither works.

The error I am getting is "Unable to open datasource 'GIS' with the following drivers", followed by a list of available drivers. Obviously it is picking up the "GIS" part from mid-way through the paths and treating it as the data source as you can see below.

The generated BAT file looks something like this:

OGR2OGR -f "ESRI Shapefile" "S:\000 GIS Files\TEST WORKSPACE 2\convert test\Converted\SHP\" "S:\000 GIS Files\TEST WORKSPACE\SiteForm.TAB" -overwrite -skiperrors


texts: Reclassifying Raster to Vectorize in ArcScan

texts:
I have a raster file featuring sidewalks for an entire city. I created this raster from polygons with the intention of using ArcScan to extract a centerline from the raster. I know that the raster has to have exactly 2 values within it. However, I'm having a problem with the vectorization where only one of the values is able to be vectorized at a time, but I cannot control it. Is the raster supposed to be in the form of, say, 0 = no data cells and 1 = sidewalks? If so, how can I change the values to match that? I've tried simply reclassifying the raster to be 0 and 1 (see below), but I get this error:

ERROR 999999: Error executing function.
Failed to execute (Reclassify).


texts: ArcGIS Flex API: Keep info window on top of other widgets?

texts:
I am working on a web mapping widget using the ArcGIS flex API and ArcGIS Flex Viewer.  I have a infoWindow popup when there is a mouse rollover on a map graphic.  Problem: the infoWindow renders below other widgets which obscures the data in the info window. Is there a way to force the info window to render on top of all the other viewer widgets?

Thanks.


texts: Projecting raster z-values from foot to meters?

texts:
I have a DEM raster (cell size 2.5*2.5). The spatial reference is lambert_Conformal_Conic_2SP and linear unit is Foot_US (0.304801).

Values for DEM was from 879.6 to 1426.78. I wanted to reproject it in UTM (cell size 10*10) and have linear unit in meters.

I simply used the project raster function and got an output raster that ranged from 879.63 to 1426.03. Does that seem right?


texts:
Let me state up front that I'm totally new to the field of GIS. I'm doing a project involving a database of buildings and their geometries in my home country (the Netherlands). At the moment, I'm trying to read geometry data from a geopackage using a python script.

If I understand the geopackage specification correctly, geometry data in a geopackage is a form of a well-known-binary. I tried to find a python package that can translate the binaries in the geopackage to usable geometries. I.e., I want to be able extract the points/linestrings/etc from the binaries.

I thought I could find a python package that would be able to do this, but I either don't understand how they work or they don't work with the data I have.

As an example, I have the following binary:

or, as represented in python:

b'\x40\x71\x00\x00\x71\x3d\x0a\xd7\x05\xc3\x0e\x41\x23\xc3\x0e\x41\x71\x3d\x0a\xd7\xe3\x1c\x22\x41\xec\x1c\x22\x41\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\xeb\x03\x00\x00\x01\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00\x71\x3d\x0a\x41\xcd\xcc\xcc\xcc\xe5\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x18\xc3\x0e\x41\x71\x3d\x0a\x41\x00\x00\x00\x00\x00\x00\x00\x00\x29\x5c\x8f\x41\x52\xb8\x1e\x85\xea\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00'

When I try to read a binary using for example the shapely package, I get an error:

from shapely import wkb
geom = b'\x40\x71\x00\x00\x71\x3d\x0a\xd7\x05\xc3\x0e\x41\x23\xc3\x0e\x41\x71\x3d\x0a\xd7\xe3\x1c\x22\x41\xec\x1c\x22\x41\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\xeb\x03\x00\x00\x01\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00\x71\x3d\x0a\x41\xcd\xcc\xcc\xcc\xe5\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x18\xc3\x0e\x41\x71\x3d\x0a\x41\x00\x00\x00\x00\x00\x00\x00\x00\x29\x5c\x8f\x41\x52\xb8\x1e\x85\xea\x1c\x22\x41\x00\x00\x00\x00\xcd\xcc\xcc\xcc\x10\xc3\x0e\x41\x5c\x8f\xc2\x41\x00\x00\x00\x00\x00\x00\x00\x00'
point = wkb.loads(geom)
print(point)

shapely.errors.WKBReadingError: Could not create geometry because of errors while reading input.

To read geometry data from geopackages in python, will I have to write my own parser?


texts: Error 193 while running gdal_calc.py using subprocess.call in Python

texts:
The following code works when used with cmd prompt:

C:\gdal_calc.py -A C:\a.tif --outfile=C:\result.tif --calc="A+1"

When I want to run the same code with subprocess:

import subprocess

subprocess.call(['C:\\gdal_calc.py', '-A', 'C:\\a.tif', '--outfile=C:\\result.tif', '--calc="A+1"'])

The following error results:

[Error 193] %1 is not a valid Win32 application

Please help me. I am highly frustating with it. This problem is basically the syntax error when the parameters are fed up into subprocess.call method what I think.

Here is the reference: http://www.gdal.org/gdal_calc.html

gdal_calc.py file is here: http://svn.osgeo.org/gdal/trunk/gdal/swig/python/scripts/gdal_calc.py


texts: Create own color ramps from attribute table based on unique column in QGIS

texts:
I have Four columns like MapCode R G B, I would like to create my own colour from the attribute table(R, G, B-Column) for Mapcode.

I would like to apply symbology based on MapCode values(which is unique) and I have the color Values (R, G, B,) as seprate Columns in the same table. how to create color ramps for unique map_code

MapCode, R, G, B

AB1C, 243, 255, 128

BB12, 124, 125, 120


texts:
In QGIS 3.4.4, I've set a field to double numeric, precision 1, to give a useful value for hectares. The field shows nn.n but the label goes to loads of decimal places.

I've tried round but it doesn't do anything (mostly). I've looked at a number of threads, including this How to change length to 2 decimal places? but it's not working.  If I set it to round only, I get the non-decimal rounded figure, but if I specify 1 or 2 places, it is ignored. I've checked & rechecked my expression [  round ("Area ha",1)  ]

Any thoughts?


texts:
I have two rasters depicting climate in two periods and I want to compare their values.

I want to create a third raster returning:

  Values of raster 1 are within the minimum and maximum values of raster 2.

and fourth raster:

  Values of raster 1 are within the 2.5 % and 97.5 % percentile of raster 2.

I have never worked with raster calculator and I do not know the syntax.


texts: arcpy.da.searchCursor - need to go through all records, but only go through one. What am I missing?

texts:
My script supposed to go through all rows in one dataset, fl1:

select one row at a time,
based on that selection, select spatially from another dataset, fl2
if spatially identical records were found, update some columns from fl1 to fl2.

My script only does through one record. What am I missing to get going through every row in fl1? I thought that SearchCursor goes from one row to the next.

with arcpy.da.SearchCursor(fc1, fields1) as search_cur:
    for search_row in search_cur:
        with arcpy.da.UpdateCursor(fc2, fields2) as upd_cur:
            for upd_row in upd_cur:
                # select this record in fc1
                whereClause = "ObjectID_1 = {0}".format(search_row[0])
                #fl1 = rcpy.MakeFeatureLayer_management(fl1, 'fl1') # creating feature layer
                arcpy.SelectLayerByAttribute_management(fl1, "NEW_SELECTION", whereClause)
                # check for spatial intersect
                #fl2 = arcpy.MakeFeatureLayer_management(fl2, 'fl2') # creating feature layer
                arcpy.SelectLayerByLocation_management(fl2, "ARE_IDENTICAL_TO", fl1, "", "NEW_SELECTION")
                # check for selection in fl2
                desc=arcpy.Describe("fl2")
                if not desc.FIDSet:
                    # if selection is not empty
                    print search_row[1]


texts: Best way to figure out which regions a point is in a shapefile?

texts:
I'm trying to figure out the best way on Android to figure out if a point in region as defined in a shapefile. For instance, I might want to know what country I'm in using the Geo Common's Shapefile file for Countries. Is there an existing Android library that can manage such tasks with Shapefiles, or an equivalent substitute (KML, for instance)? A few other points:

I might be using a lot of regions. For instance, I might want to know what county I'm in, along with state and country. Thus, there is potential for overlap of regions.
Speed is desirable, if possible. I might, for instance, figure out which country I'm in, then state, then county, and only look at the applicable states for a country, and counties for a state.
Compactness is very important. At least one of the data files I want to work with is 47 Mb as a KML, Shapefiles seems like a reasonable size, assuming I can find a decent library for it.


texts: Polygons get stretched all over the map when changing CRS in QGIS

texts:
I am trying to change the CRS of a vector layer in QGIS, but the output seems to be deformed.

As you can see from the attached pictures, some of the points belonging to polygons which are situated near the edges of the map in the input layer get wrapped on the other side of the world in the reprojected layer. I also wish to point out that I'm properly using "save as...", and not "set layer CRS".

Can anybody explain why is this happening? Is it because, as I suspect, some points are exactly on 180 longitude and not, for instance, 179.999?

Is there is a way to correct the problem (possibly without manually editing each problematic point)?


texts: Reproject image collection according to the projection of another image collection

texts:
I have two image colelctions. one of sentinel 1 and one of sentinel 2.
I want to reproject the sentinel 1 so its' projection and cell size would be like the sentunel 2 images.
It worked for me when I did it without function on only one image, but now when I have ttried to create a function for this, when I check the results I see that the pixels of the two polygons do not sit on top of each other, e.g seems like the reprojection didn't work.

This is my function:

////Fix reproject SAR 1////

var SARresample=Clipped_Speckles.map(function(im){
  var sar1Projection =  im.projection();
  var NDVIProjection =  ee.Image(listOfImages.get(listOfNumbers[0])).projection();
  var SARreproject=im.reduceResolution({reducer: ee.Reducer.mean()}).reproject({crs: NDVIProjection});

   return im;
});

after the function. I have tried to display two images in order to see the pixels but it seems like the function didn't function:

//Visualize the resample vs NDVI

var listOfIResamples =(SARresample.toList(SARresample.size()));

var NumberOfImagesSAR=[0];

var ResampledIM= ee.Image(listOfIResamples.get(NumberOfImagesSAR[0]));

Map.addLayer(ResampledIM,{bands: 'VH_Filtered',min: -30, max: 1},'SAR');

//visualize NDVI

var listOfImagesNDVI =(filterNDVI.toList(filterNDVI.size()));

var NumberOfImagesNDVI=[0,1,2];

var ndviIM= ee.Image(listOfImagesNDVI.get(NumberOfImagesNDVI[1]));

Map.addLayer(ndviIM, NDVIcolor,'NDVI');

Here you can see the full code:

https://code.earthengine.google.com/0b322b779162bc1eb18a6143867d1a99

*I have changed the speckle filter size to 2 for it to be easier to see the pixels.

*the polygon is in north of Thailand , cannot share it but any random polygon in that area will work.


texts: Custom Geoprocessing messages/progress indicator Arcgis 10.2 Modelbuilder?

texts: I was wondering if anyone knew a way of customizing the progress indicator for tool developed in Arcgis's model builder. I have a complicated tool, which uses a feature selection iterator, and I would like to either indicate what percent of features have been processed, or output a message with the current count.  I haven't found anything for model builder, but arcpy has the AddMessage tool.


texts: How to check integrity of Python-processed shapefile

texts:
In my QGIS project of American railroads, using USGS Base Topo Layer and shapefiles giving railroad details as of 2015, I have used the Python-3 pyshp library to edit the dbf database, adding some fields.
After writing out the edited shapefile with pyshp.Writer QGIS rejects my new shapefile as invalid, without any pointer as to the reason of rejection.
How can I check the integrity of a processed shapefile, preferably with Python ?
I use QGIS-3.14.11-Madeira, pyshp-1.1.12 and Python-3.6


texts: producing smooth iso polygons from raster surface?

texts:
The objective is to turn a raster surface (double value type) into smooth, visually appealing polygons representing iso-levels. I cannot stress enough that these must be visually attractive. Jaggedy lines that are raster remnants are not acceptable. Here are the parameters of my problem:

Input: Raster surface of atmospheric forecasts (not a DEM). Z values can fluctuate greatly over a short horizontal distance, but the surface does tend to have sane looking clusters of pixels with peaks and valleys.
Desired output: Smooth polygons representing iso-levels at irregular intervals. Class breaks are 0-5, 5-20, 20-75, 75-150, and >150.
I really have artistic license to make these look nice. I have been requested to remove all remnants of jagged raster edges. Obviously, I want to be as true as possible to the underlying data, but the marching orders are something that is simple and easy for the novice user to understand without a lot of details distracting the eye.
Needs to be automated to run on regular intervals. I can code in several languages, but prefer python or command line utilities.
I have access to ArcGIS or can use FOSS tools, but prefer to run in Linux, which probably favors the latter.

I'll talk a bit about what I've tried below, but everything I try runs into problem after problem, to the point where it's just getting too complicated. My inner mental dialog is saying "It's 2017, somebody must have really solved this problem by now!" So... am I missing the obvious solution?

First, I tried reclassifying the raster into categories, and vectorizing that to polys. I got similar results from Esri "reclass by ASCII file" and gdal_calc. The big problem here is that it produced lots of areas where there adjacent polygons skipped a category or meet in diagonal corners (I've flagged a few in the diagram), and I'd really like concentric bands. I might be able to reduce the cell size and smooth this enough to solve the problem, but I'm not sure how far down this path I want to go...

Even with the polygon simplification in ArcGIS, I still had to smooth the lines to make them look halfway decent, but this produced some strange topology issues with holes in my coverage. I'm sure I could solve this, too, but I feel like I'm just going down multiple rabbit trails at this point.

The other thing I tried was to create contours. Esri's "contour list" tool produced slightly nicer results than the gdal_contour command, but I was able to produce contour lines using both. Ultimately, those were just lines, not polygons. I found a few references to converting to polys, but this feels like I'm re-inventing the wheel, so I wanted to come up for air and see if this is on the right track. (Esri toolbox produced error when run, blog post not automated yet)

There is also the Esri "GA Layer to Contour" tool, but I was having trouble figuring out the inputs to that. It was unclear to me how to create the input layer needed for that tool and the documentation was less than clear and lacking useful examples. It looks like I need to convert my raster to points and interpolate it again? Any hints on this one?

I feel like I have a number of really complicated options, but this seems like a problem that should have been solved a long time ago. I also feel like I have a lot of useful tools, but I'm not putting them together in the right order...


texts:
I'm a (relatively) new Mac QGIS user, and I'm having trouble with my processing toolbox.

Specifically, I'm having an issue with the Orfeo toolbox not activating after downloading the last update (2.18.14). I successfully re-downloaded the Orfeo Toolbox folder from this website: https://www.orfeo-toolbox.org/download/

It downloaded into one of the "users" on my computer, but I moved the file to the Applications section where QGIS downloaded. When I open Processing>>Options I see Orfeo toolbox, and it's checked as being active. See below:

I've checked that the applications and command line files are correct from other message boards on this website, so I don't think that should be the issue at this point but let me know if they're not correct.

Under the menu section of the Processing>>Options section, I see the Orfeo toolbox but there aren't any dropdown options like the other menu options that are actually showing up.

Any ideas as to what I've done wrong here?


texts: How to get QGIS to display grayscale of an SRTM geotiff?

texts:
I cannot get QGIS to recognise my SRTM geotiff as having any grayscale elevation data. In the screenshots below, opening the geotiff file in MicroDEM by default hillshades and pseudocolors the hypsometry, so the file definitely holds z-value data.

The screenshot from QGIS further below, where I have opened the same geotiff file and checked the properties histogram, shows that none of this elevation data is apparently recognised by QGIS. I must be doing something wrong, but what?

(SRTM = Shuttle Radar Topography Mission)

My geotiff opened in MicroDEM (good):

The same geotiff opened in QGIS (bad):


texts: How to connect points to lines in order using QGIS?

texts:
I have a set of points and I want to create lines from them. Points belonging to one line have unique UniqueID to be grouped by. I try to use QGIS' Points to path algorithm. What I don't have is an order field, so the problem is that it connects the points in random order.

The picture contains an example of a line with 6 points. The QGIS Points to path algorithm connects the points in the following order: 1,3,2,6,4,5 and creates visually a good looking line, but with the length longer than the line between point 1 and 6. So I would like to know how to establish the order for the points to be added when creating a line. The points need to be connected from beginning to end in order i.e. 1,2,3,4,5,6.

The attribute table is shown on the picture. Each point has unique fid and a field called UniqueID that groups the points to one line i.e. points that belong to the same line have the same UniqueID. The name of the attribute comes from the data I have received so can't really make it better.
How can I introduce an order so that I can have points connected by shortest path? The end goal is to calculate the length of lines, so if that could be done easily on the lines created from the unordered points, it is also a good solution to the issue.


texts:
While there are functions to convert geography to GeoJSON:
text ST_AsGeoJSON(geography geog, integer maxdecimaldigits=15, integer options=0);
I do not see any way to go the opposite direction (GeoJSON to geography) only to geometry.
http://postgis.net/docs/reference.html#Geometry_Constructors
Is this on the roadmap?

The data in the database has geography polygons or multi-polygons to describe each item in the DB. The query is coming from lealfet.js, which allows me to easily get the shape in GeoJSON.
What I am doing is just a simple ST_Intersects() to check for overlap and return. Yes, I have a hack to convert this to EWKT right now, but I was wondering if it was on the roadmap to support this in a future version of PostGIS, or where one adds it as a wish list (since symmetry is a good thing, and it seems a shame that only geometry is supported in that direction).


texts: Cannot load tiles from Mapbox using leaflet.js

texts:
I am a new user and I can't seem to get the tiles on the map.

var mymap = L.map('mapid').setView([51.512, -0.104], 1);

L.tileLayer('https://api.mapbox.com/styles/v1/mapbox/{mapbox.mapbox-streets-v8}/tiles/{z}/{x}/{y}?access_token={accessToken}', {
    attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © Mapbox',
    maxZoom: 18,
    id: 'mapbox.mapbox-streets-v8',
    accessToken: 'accessToken.here'
}).addTo(mymap);

var source = new EventSource('/topic/twitterdata2');

source.addEventListener('message', function(e){
    obj = JSON.parse(e.data);
    console.log(obj);
    lat = obj.place.bounding_box.coordinates[0][0][1];
    long = obj.place.bounding_box.coordinates[0][0][0];
    username = obj.user.name;
    tweet = obj.text;

    marker = L.marker([lat,long],).addTo(mymap).bindPopup('Username: ' + username + '
Tweet: ' + tweet + '');

}, false);


texts:
I am using ArcGIS 10 and I would like to select (in a manner of well spatially distributed) some of the polygons included in a shapefile as testing population for a statistical analysis.

The rest of the polygons included in the shapefile I intend to use them for the validation of the statistical outcome. As I am not so familiar with statistics, is there a way to automatically select a well spatially distributed sample in ArcGIS or if not, can you please suggest other solutions?


texts: Drawing buildings with Holes inside with Openstreetmap XML data

texts:
I'm trying to draw a map using OSM raw data.Some of buildings have holes inside (may even more than one hole at a time) .

In way list, all I can see is set of "way nodes". Where should I get the information in OSM.xml that a  Node represent a hole is belong to a specific building?

Note: In Relation there is innerWay - outerWay relation for some buildings but it does not cover all the cases as I see. How openstreetmap draw building?

I achieved drawing it in Opengl using tesselation without knowing which hole belong which building but, I couldn't apply same approach on unity. Here is a screenshot from opengl :


texts:
I am trying to correct ellipsoid heights to orthometric height, transforming from WGS84 to EGM2008 using the EPSG code 3855. I have tried with pyroj, as in this question:
import pyproj
EGM2008height = pyproj.Proj("EPSG:3855")
wgs84_ll = pyproj.Proj("EPSG:4326")
print(pyproj.transform(wgs84_ll, EGM2008height, 34.68016909181368, 38.31245226053967,100))

which returns the error:
(Internal Proj Error: proj_create: unrecognized format / unknown name)

I understand that I need the egm08_25.gtx file, which I have downloaded but I don't know how to use it so with pyproj. Is there a way to add EPSG codes to pyproj?
I can get it to work by looping over my data and using gdaltransform by calling the command:
cmd = 'echo ' + '"' + str(X[i]) + ' ' + str(Y[i]) + ' ' + str(Z[i]) + '"' + ' | gdaltransform -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/Usegm08_25.gtx"'

But this is a slow process and I think using pyproj will be more efficient.
There are several questions on this already, but none of them comprehensively show how to use the .gtx file or add EPSG codes into pyproj.
EDIT.....
Related Qs where I think the answe requires reviewing:
Vertical Datum transformation using Pyproj
Converting Ellipsoidal height to Orthometric height in Python
Converting EPSG:2263 to WGS84 using Python Pyproj


texts: Spatially clipping Polygons using ArcMap

texts:
I am using ArcMap.

How can I spatially clip the polygon shapes outlined in red marking so that the blue polygons outside the red perimeter are deleted (red marking is just made with snipping tool pen)?

The clipped blue polygons need to share the same boundary as the purple polygons and there are gaps between some of the "islands" of purple polygons (i.e. they are not connected). Also, the blue and purple polygons are in 2 different layers (shapefiles) and there is overlap between them.


texts: Looking for bathymetry data for Lake Sakakawea in North Dakota

texts: Navionics has a nice web map that has some incomplete coverage of Lake Sakakawea. Garmin's BlueChart looks to have something for about $160, but appears to be only compatible with GPS units. Don't need the entire lake just a few areas of it. Would like the data to be ArcGIS compatible. I have attempted to do raster to vector in ArcGIS taking screen shots of Navionics web map but haven't had much success. Any ideas / help would be most appreciated.


texts: Is it possible to set up a local QGIS Plugin repository?

texts:
I've developed a couple of custom QGIS Plugins for use by our office. I can't post them to the official QGIS Plugins repository as they reference/utilize sensitive data on our own servers.

For distribution, I'm currently just having users copy the plugin folder over to their local plugins folder. However, this is inconvenient and won't give them an easy way to update.

Can I set up a custom repository on a network drive or server (behind our firewall) that users can connect to?


texts: Multipolygon in Los Angeles County Parcels Tax Roll

texts:
I downloaded the Los Angeles Parcels Tax Roll data into QGIS. I copied the information for a single parcel results in QGIS:

wkt_geom   OBJECTID_1  AIN Roll_Year   AssessorID  TaxRateArea TaxRateArea_CITY    SitusHouseNo    SitusFraction   SitusDirection  SitusStreet SitusUnit   SitusConcatenated   SitusCity   SitusZIP    SitusZIP5   UseCode UseCodeDescChar1    UseCodeDescChar2    UseCodeDescChar3    UseCodeDescChar4    totBuildingDataLines    YearBuilt   EffectiveYearBuilt  SQFTmain    Bedrooms    Bathrooms   Units   RecordingDate   Roll_LandValue  Roll_LandBaseYear   Roll_ImpValue   Roll_ImpBaseYear    Roll_totLandImp Roll_PersPropValue  Roll_FixtureValue   Roll_HomeOwnersExemp    ExemptionClaimKey   ExemptionType   ExemptionFullOrPartial  Roll_RealEstateExemp    Roll_PersPropExemp  Roll_FixtureExemp   netTaxableValue__Land_Imp_REX_HOX_  totBoundaryDescLines    BoundaryDescLine1   BoundaryDescLine2   BoundaryDescLine3   BoundaryDescLine4   BoundaryDescLine5   BoundaryDescLastLine    ParcelClassification    AdministrativeRegion    CENTER_LAT  CENTER_LON  Shape_Length    Shape_Area
MultiPolygon (((
6445171.31999936699867249 1754319.45000045001506805
6445165.38999937474727631 1754259.7400003969669342
6445040.50999936461448669 1754272.15000040829181671
6445043.06999936699867249 1754297.81000044941902161
6445046.4499993622303009 1754331.86000040173530579
6445171.31999936699867249 1754319.45000045001506805
))) 539190  7514021016  2014    7514-021-016    9359    TORRANCE    249         PASEO DE LAS DELICIAS       249 PASEO DE LAS DELICIAS   TORRANCE CA 90277-6439  90277   0100    Residential Single Family Residence Unused or Unknown Code (No Meaning)     1   1948    1948    1456    3   2   1   19900220    70378   1975    34771   1975    105149  0   0   7000                0   0   0   98149   1   TRACT # 10302 LOT   16 BLK    N                         14  33.81230373 -118.38371317   370.98874200978423  7529.721003203928

What are the units of the multipolygon? My QGIS is set to EPSG: 3857. Although I think the units for the MultiPolygon are the same no matter what CRS is used.

And then how do I get just this parcel into QGIS. The database is rather large and I don't need it all.


texts: Dealing with NoData in ArcGIS Spatial Analyst?

texts:
When using Spatial Analyst I regularly have issues with NoData values cropping up in the data.

Various tools appear to introduce NoData values regularly.

These values then cause issues with other tools and have to be eliminated with reclassification, usually resulting in other problems.

How should NoData values be dealt with?  Is there a simple way to replace NoData values with zero?  Can NoData be avoided altogether?


texts: Georeferenced image not displaying in QGIS

texts:
I'm still on my 'L' plates. I'm trying to georeference a map in QGIS so that I can trace out several features. I've done this (once) before, and don't understand why it wont work now.
I've a pdf map, that I've converted to tiff and jpg. They all load into the georeferencer. I've got 20 or so good control points spread across the map
I've set the projection to the same as my training file. The georeferencer will run, creates a new image file and loads this into QGIS. QGIS will not display the image.
The image is legit, and graphics software will display it.
QGIS will not display the new layer - even though the layer is loaded.
'Zoom to layer extent' shows a blank, white display.
Does anyone have any suggestions?


texts: Axis of approximate symmetry in irregular polygons

texts: Does anyone have an idea how to determine the axis of approximate symmetry in irregular convex polygon? (Axis of approximate symmetry being defined as straight line cutting the polygon in two halves in a way that their shapes are as similar as possible.)


texts:
I am trying to add a jpeg200 image to my postgresql database.
Problem: this format is not implemented by default in postgresql/postgis which uses the gdal library.

After some research, I found that jp2 needs another library http://www.gdal.org/formats_list.html (named libjasper). I've managed to get a dll of it but I have no clue where to put it, so far I've tried in \bin and \lib in the postrgresql's repository but no luck so far.

Then, I've found I need the dll compiled with the same compiler as postgis.

Has anyone been confronted with this problem before?

I could make my way around by converting my jp2 into tiff (ogr, qgis, arcgis) but as you may know tiff are way bigger than jpeg2000.


texts: PyQGIS processing: Overwrite not possible

texts:
I have a QGIS plugin that runs the v.clean algorithm two times (first time snap, second break) on a chosen layer:

    output = self.dlg.directory.text()

    params_snap = {
        'input' : input,
        'type' : 1,
        'tool' : 1,
        'threshold' : 0.000095,
        'output': output,
        'error' : error,
        '--overwrite': True,
        'GRASS_OUTPUT_TYPE_PARAMETER': 2
    }
    processing.run('grass7:v.clean', params_snap)

    params_break = {
        'input': output,
        'type': 1,
        'tool': 0,
        'threshold': 0.1,
        '-c': True,
        'output': output,
        'error': error,
        '--overwrite': True,
        'GRASS_OUTPUT_TYPE_PARAMETER': 2
    }
    processing.runAndLoadResults('grass7:v.clean', params_break)

My problem is that after running snap, the created output file is not overwritten. runAndLoadResults of the break algorithm loads the output file as expected, but without broken lines. So only snap was applied. Since --overwrite is True it should work, or do I miss something?

If I choose another output file for the break algorithm everything runs perfect, so the problem is not the algorithm itself.

As stated in my other question here, a temporary output is not possible for Grass algorithms and therefore no solution.

EDIT: The processing log messages print an error:

2019-01-30T11:59:16     INFO    C:\Windows\System32>v.out.ogr -c type="line"
input="output6c3d7544cb524a2486d97aca71a41bc4" output="C:\Users\USER
\Desktop\test.gpkg" format="GPKG" --overwrite
2019-01-30T11:59:16     CRITICAL    ERROR 1: A file system object called
'C:\Users\USER\Desktop\test.gpkg' already exists.

So Grass recognizes overwrite, but it does not work. Is this the wrong command, for what I try to achieve? The algorithm help (processing.algorithmHelp("grass7:v.clean")) gives no information about overwrites. I got the information from the v.clean manual.


texts: Longitude/Latitude from NetCDF with (row,column)

texts:
I downloaded a "Sentinel-3 Land Surface Temperature" netcdf file. I want to get the longitude and latitude for each temperature entry in the file. Well the file has pairs of rows and column. Using Python I can get the temperature using :

from netCDF4 import Dataset

with dataset as ncFile:
    nc_attrValue = ncFile.variables['LST'][1000][1000]
    print(type(nc_attrValue))

So any idea how I can convert [1000][1000] to longitude or latitude. Here is a snapshot of the file netcdf file converted to .cdl using NCO:

    netcdf LST_in {
dimensions:
    columns = 1500 ;
    orphan_pixels = 187 ;
    rows = 1200 ;
variables:
short LST(rows, columns) ;
        LST:_FillValue = -32768s ;
        LST:add_offset = 290.f ;
        LST:long_name = "Gridded Land Surface Temperature" ;
        LST:scale_factor = 0.002f ;
        LST:standard_name = "surface_temperature" ;
        LST:units = "K" ;
        LST:valid_max = 32767s ;
        LST:valid_min = -32767s ;
// global attributes:
        :absolute_orbit_number = 12616U ;
        :comment = " " ;
        :contact = "eosupport@copernicus.esa.int" ;
        :creation_time = "20180720T144045Z" ;
        :history = "  2018-07-20T14:40:45Z: PUGCoreProcessor /data/ipf-s3/workdir47/568114279/JobOrder.568114279.xml" ;
        :institution = "SVL" ;
        :netCDF_version = "4.2 of Mar 13 2018 10:14:33 $" ;
        :product_name = "S3A_SL_2_LST____20180720T124923_20180720T125223_20180720T144045_0179_033_323_3240_SVL_O_NR_003.SEN3" ;
        :references = "S3IPF PDS 005 - i2r7 - Product Data Format Specification - SLSTR, S3IPF PDS 002 - i1r6 - Product Data Format Specification - Product Structures, S3IPF DPM 007 - i1r8 - Detailed Processing Model - SLSTR Level 2" ;
        :resolution = "[ 1000 1000 ]" ;
        :source = "IPF-SL-2 06.13" ;
        :start_offset = 21649 ;
        :start_time = "2018-07-20T12:49:23.587859Z" ;
        :stop_time = "2018-07-20T12:52:23.279310Z" ;
        :title = "SLSTR Level 2 Product, Land Surface Temperature measurement" ;
        :track_offset = 998 ;
data:

 LST =
  _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, -7188, -6706, -6884, -6249, -6173, -6721, -6924,


texts: Geometry checker not removing overlap between features

texts:
I have a polygon shapefile layer (CRS EPSG:27700 British National Grid) that has a deliberate overlap and gap between features.

When using the Geometry Checker Plugin (V0.1) in QGIS 3.10.9 I am selecting the topology checks 'check for overlaps' and 'check for gaps' both with 100 map units sqr.

The geometry checker finds both errors.

I select to fix both errors with the default resolution and the geometry checker goes through the process of fixing the errors, telling me they are fixed and generating a new shapefile.

However, the overlap between two features in the polygon still remains.

Why isn't this overlap being removed?


texts: Calculating total area of pie charts in ArcMap

texts:
Background: I am scaling the size of pie charts by a field which indicates a percentage. I am doing this for two datasets displayed on two maps and would like the resulting areas of the scaled pies to be equal between the two maps. I think that I would be able to do this by normalizing the scaled pie sizes using one of the two total areas of the pie charts.
Question: Is it possible to find out the total area of the pies as they are displayed on the map? I am able to scale the sizes manually for example, by setting a larger minimum size symbol, however it is difficult to visually gauge equal total areas this way.


texts:
I am working in a attribute table trying to format my numbers, but want to keep ending zeros. I am using format_numbers("column",sig figs). For the most part, it is working, but for example row two in the picture below, I need to zero to show to have everything line up perfectly. I am working in QGIS 3.10.4.
Any suggestions on how to do this?


texts: Issue on Uploading Two SOE to ArcGIS Server

texts:
I need to upload two REST SOE to ArcGIS Server. They are technically same and I just copied second one from the first and renamed the project and built it but when I tried to upload the second SOE ArcGIS server Complained that

so I tried to change some configurations at second SOE as

1- in Properties/AssemblyInfo.cs I changed

[assembly: AssemblyTitle("GeometricNetworkUtility")]
[assembly: AssemblyDescription("Geometric Network Utility")]

[assembly: AddInPackage("ArcGIS.Soe.Rest.GeometricNetworkUtilitybyPhase", "f15446f6-4217-4c50-a970-078c0b64a20a",
    Author = " ",
    Company = " ",
    Description = "Geometric Network Utility",
    TargetProduct = "Server",
    TargetVersion = "10.4",

    Version = "1.3")]

To

[assembly: AssemblyTitle("GeometricNetworkUtility2")]
[assembly: AssemblyDescription("Geometric Network Utility 2")]

[assembly: AddInPackage("ArcGIS.Soe.Rest.GeometricNetworkUtility2", "f15446f6-4217-4c50-a970-078c0b64a20a",
    Author = " ",
    Company = " ",
    Description = "Geometric Network Utility 2",
    TargetProduct = "Server",
    TargetVersion = "10.4",

    Version = "1.3")]

and in GeometricNetworkUtility.cs  I changed the

Description = "Geometric Network Utility",
DisplayName = "Geometric Network Utility",

to

Description = "Geometric Network Utility 2",
DisplayName = "Geometric Network Utility 2",

but still after rebuilding the solution and uploading the new SOE I am getting the same error message. Can you please let me know how to solve this?


texts: Haversine formula on rotated poles sphere

texts:
I'm working with climate data from the CORDEX domain WAS-22. The data are on rotated poles grid defined as:
char rotated_pole;
  :long_name = "coordinates of the rotated North Pole";
  :grid_mapping_name = "rotated_latitude_longitude";
  :grid_north_pole_latitude = 79.94999694824219; // double
  :grid_north_pole_longitude = -123.33999633789062; // double

I'm trying to calculate the relative vorticity using wind components and since I don't want to lose the natural variability of the simulation output, the interpolation/projection to the lonlat grid is not an option. So I'm using the finite difference approximation to calculate the relative vorticity as:

For that, I need to calculate the delta x and delta y on the sphere of the rotated poles. I know that Haversine formula works in case of the lonlat grid and accounts for the longitudinal convergence in the direction of the poles. My question is how the Haversine formula can be modified to work in the case of rotated poles or if there is another option?


texts: Adding delimited text layer with Python outside of QGIS

texts:
I'm trying to create a QGIS project using Jupyter environment.
Until now I managed to add a Raster layer.
I now want to add a csv file, in QGIS I would add it like that: Layer -> Add Layer -> Add delimited Text Layer.
How can I do the same using Python?
This is what I have now:
from qgis.core import QgsProject,QgsRasterLayer,QgsRaster,QgsVectorLayer
import boto3
from PyQt5.QtCore import QFileInfo
from qgis.gui import QgsMapCanvas,QtCore
from PyQt5.QtWidgets import QApplication, QWidget,QMessageBox,QDialog

project = QgsProject.instance()
#here I add the raster layer
fileInfo = QFileInfo('rasterlayer')
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer('rasterlayer', basename)
if not rlayer.isValid():
   print("Layer failed to load!")
QgsProject.instance().addMapLayer(rlayer)

This is what I tried to use to add the delimited text layer and I keep geeing Invalid error.
uri='file:///coordinatestry.csv?delimiter=,&yField=y&xField=x'
layer = QgsRasterLayer(uri, 'name', 'delimitedtext')

uri = ("file:///{}?"
               "&delimiter={}"
               "&xField={}"
               "&yField={}"
               "&crs={}"
               "&decimal").format('coordinatestry.csv',
                                        ",",
                                        "x",
                                        "y",
                                        "EPSG:4326",",")
vlayer = QgsRasterLayer(uri, "layer", "delimitedtext")

fileName = 'coordinatestry.csv' # the path to the csv
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
uri = "file:/"+fileName+"?delimiter=%s" % (",") # tried without the file:/ too and different delimiters
vlayer = QgsVectorLayer(uri, baseName, "delimitedtext")

I keep getting this: <QgsRasterLayer: 'layer' (Invalid)>
The end goal is to use Atlas export as images for the coordinates in the csv on the raster layer.
Any ideas how to make the adding delimited text layer work?
Edit: Tried this as well as was suggested in comments:
uri = "file://{}/coordinatestry.csv?delimiter={}&xField={}&yField={}".format(os.getcwd(), ",", "inside_x", "inside_y")
vlayer = QgsVectorLayer(uri, "layer name you like", "delimitedtext")

Still got : <QgsVectorLayer: 'layer name you like' (Invalid)>


texts: QGIS 3.6: Vector "Polygons to Lines" missing

texts: I use the tool Vector-->Geometry Tools-->Polygons to Lines all the time.  It is in the menu up until 3.4. In 3.6, that menu item is no longer available.  Where did it go?  Is there another way to do this in 3.6?


texts:
Suppose I have this list of rasters:

(L <- setNames(replicate(3, raster(matrix(runif(100), 10))), c('A', 'B', 'C')))

How can I unlist a list of rasters to individual objects in R?

With a list of dataframes, I can use something similar to the following code:

df <- data.frame(matrix(unlist(l), nrow=132, byrow=T),stringsAsFactors=FALSE)


texts: Converting large database of addresses and associated data into point shapefile within ArcMap

texts:
I have a large Excel spreadsheet of addresses and project data. I would like to automatically convert this into a point shapefile based on the addresses.

Since these are not coordinates and are addresses, how do I do this using ArcMap?


texts: Calculating mean value of 280 fields with Field Calculator from attribute table of shapefile?

texts:
I have to average 280 fields in 1 field. Each field has 300 rows.

That is my problem. I think I need a python script. I don't want to create a raster for each field and then use the raster calculator in ModelBuilder with iterations...

I don't know how to make that script in Python.

In this pic, the field "Prom_P2" is what I need to fill with the mean or average of the previous 280 fields.

This is code that @BERA shared with me, and it works:

import arcpy
shapefile = r'D:\PYTHON_practica\pp_1998_2015.shp'

field_names=[f.name for f in  arcpy.ListFields(shapefile,wild_card='*Lay*')] #Create a list of all the fields with 'Lay' in field name
field_names.append('promedio') #add your mean field to the list

with arcpy.da.UpdateCursor(shapefile,field_names) as cursor:
 for row in cursor:
     row[-1]=sum(row[:-2])/len(row[:-2]) #Calculate mean and assign this to the mean field (-1)
     cursor.updateRow(row)


texts: How to use offline routing Dijkstra in javascript?

texts:
I am using a Openlayers 3 to create a Hybrid mapping app and among the constraints it must be for offline use, I have a road network in a geojson file and I need to get the shortest path from the location of the user to a destination following the roads network.

I tried pgrouting before and it worked well but now I need to do so without relying on any server.

Is there any way I can make it with javascript and Openlayers 3? Is there any library that can do offline routing ?


texts: Clip in QGIS 2.4.0 resulting in empty output files

texts:
I'm seeing the same problem to Clip or Intersection tool - QGIS giving empty files trying to use the Clip tool since upgrading to QGIS 2.4.0, but have not been able to fix the problem based on the responses there. I tried creating a simple example to ensure was not my data but am seeing the same problem.

I have a series of points & a polygon.  I want to retain those points within the polygon only.

I use Vector -> Geoprocessing Tools -> Clip

This results in an empty attribute table.

Options I've tried:

I have tried switching the Input and Clip layer in case I had these round the wrong way
Disable on the fly reprojections (confirmed the layers do overlap, they are in the same projection)
Checked through Properties to confirm both layers have CRS EPSG:4326 - WGS 84
I left this simple example running for approximately 40 minutes with no change in outcome in case I was just not leaving enough time (in previous versions this was successful in a short time span).
I have tried using other Geoprocessing tools to see if it is just clip - I do not receive any output with any of them.


texts:
I'm learning how to use Oblique Mercator projection in QGIS 3.4.7. I created a custom projection with "omerc" and passed two points (New York City, USA and Tokyo, Japan) as the parameter.

Extent: Extent not known
Proj4: +proj=omerc +lat_1=35.652832 +lon_1=139.839478 +lat_2=40.728157 +lon_2=-74.077644 +ellps=GRS80

I loaded the country shape file from Natural Earth to layer.

However, when I apply my custom projection, the map has distortion at the edges. The edges of the maps are unexpectedly connected.

Can somebody share how to solve this problem?

I found that proj4 omerc demo image also has the same problem of zingers.


texts: How to get Subtype name in Field using field calculator

texts:
I would like to get feature subtype name in another field using field calculator.

Any one can advise to make this success?

I am expecting ".Cal" file code using with Field calculator.


texts: How can I insert a record N amount of times using cursors?

texts:
I need to generate a polygon FC using some records previously flagged as 'NUMBER_GT_NUMBER_END = 1'

These records need to be loaded into a new schema of the same FC in order to include duplicates that will be used later for other processes. If you run the code as I am giving it here, it will print out the number of records per OBJECTID. The issue is that I am able to insert only one for each case. Once added a record into the new FC schema, I want to update the SITUS_STREET_NUMBER field with the value stored in x + 1. I have been going back and forth with other approaches like using Append, Copy, among others, but the process is really slow in comparison to the one achieved so far using the da.cursors.

Example:

OBJECTID = 28964 ---> 3 records

4815

4816

4817

This means that the above polygon geometry with OBJECTID = 28964 needs to be inserted to the new FC three times and that the SITUS_STREET_NUMBER field needs to be updated with the listed values.

Script:

import arcpy, re

arcpy.env.workspace = r'YOUR DIRECTOTY PATH HERE\locator_address_ranges.gdb'
arcpy.env.overwriteOutput = True

# Set a variable for the workspace
workspace = arcpy.env.workspace

featureClass = workspace + '\\' + 'ParcelsTax_Formated'
featureClassCopy = workspace + '\\' + 'ParcelsTax_Formated_new'
number_Field = "SITUS_STREET_NUMBER"
number_end_Field = "SITUS_STREET_NUMBER_END"
oid = "OBJECTID"
number_Field2 = "SITUS_STREET_NUMBER2"  # Long Integer
number_end_Field2 = "SITUS_STREET_NUMBER_END2"  # Long Integer
countAddrPnts = "Points"  # Long Integer
full_addr = "SITUS_FULL_ADDRESS"
greater_than = "NUMBER_GT_NUMBER_END"  # Value = 1 --> SITUS_STREET_NUMBER2 is > than SITUS_STREET_NUMBER_END2 - Short Integer
fieldsToUpdate = (number_Field, full_addr)

dsc = arcpy.Describe(featureClass)
fields = dsc.fields

# List all field names except the OID field

fieldnames = [field.name for field in fields]
fieldnames.append("SHAPE@")
print "Updated List : ", fieldnames

## STEP 1
# Insert rows for scenario where SITUS_STREET_NUMBER2 is > than SITUS_STREET_NUMBER_END2
arcpy.MakeFeatureLayer_management(featureClass, "countAddrPnts_lyr")
arcpy.SelectLayerByAttribute_management("countAddrPnts_lyr", "NEW_SELECTION", 'NUMBER_GT_NUMBER_END = 1')

rows = arcpy.SearchCursor(featureClass, 'NUMBER_GT_NUMBER_END = 1')
row = rows.next()

while row:
    arcpy.SelectLayerByAttribute_management("countAddrPnts_lyr", "NEW_SELECTION",
                                            "OBJECTID = " + str(row.getValue(oid)))
    print "OBJECTID = " + str(row.getValue(oid)) + ' ---> ' + str(row.getValue(countAddrPnts) + 1) + ' records'

    with arcpy.da.SearchCursor(featureClass, fieldnames, 'NUMBER_GT_NUMBER_END = 1 AND OBJECTID = ' + str(row.getValue(oid))) as sCur:
            for x in xrange(int(row.getValue(number_end_Field2) - 1), int(row.getValue(number_Field2))):
                print x + 1
                with arcpy.da.InsertCursor(featureClassCopy, fieldnames) as iCur:
                    for row in sCur:
                        iCur.insertRow(row)
    row = rows.next()

GIS data (FileGDB) and python script here


texts: Using wildcard to Append multiple feature classes in ModelBuilder?

texts:
I have a folder of shapefiles and each one starts with a code.  There are multiple shapefiles that start with each code, and I would like to create a model in ModelBuilder that will select the shapefiles that have the same code, and append them into the appropriate existing geodatabase feature class.

Is there a way to use a wildcard in the Append tool, like to specify that any and all files in a certain folder starting with 'XPGC*' get appended into the 'Pipe_Grade_Change' feature class?  Or is there a way to use an iterator to accomplish this?


texts: Change size of frame for qgis map tip containing text

texts:
I have been trying to increase the size of my maptip frames on MacOS using Safari with limited luck. I have changed the size of the frames but have not been able to display text inside the frames.  My knowledge of html is pretty basic .. about nil.

I have tried:

[% "LIGHT_COND"  ||  '\n'  || 'Persons ' ||  "TOTAL_PERS" %]
<iframe width=3000 height = 3000>
</iframe>


texts: Changing field values with ArcPy update cursor?

texts:
I want to figure out how to be able to change values in fields for use in further projects. Right now , I would like to be able to add '-WA' to the existing field value 53. (ie 53-WA). This is a test project as I am trying to teach myself python. The code I got so far runs but the values do not change.

import arcpy, os

arcpy.env.workspace = r'C:\Users\david.fleck\Python.mdb'

fc= 'Parking_Lots_Collected_2015'

field = ['STATEFP']

cursor = arcpy.UpdateCursor(fc,field)

with arcpy.da.UpdateCursor(fc,field)as cursor:

    for row in cursor:
        if row== 53:
            row.setValue(row.getValue(field) + '-WA')
            cursor.updateRow(row)

print "Finished"

I have also tried row ==''53 + '-CA' instead of row.setvalue but still no change.


texts: Create new vector of polygons that are within 500m of a point in vector of points P

texts:
If I have a vector of polygons A (say 200) and want to create a new vector (say 120) of the original 200 polygons that are with 500m of a vector of points P (say 50), this is the method I do. I want to know if there is a simpler way

1.) Add calculated field to A called "rownum" with expression $rownum to create an ID field for the polygons

2.) Create buffer vector B of original 200 polygons from A with a 500m radius

3.) Do a spatial query which does

   source     : B
   operator   : INTERSECTS
   references : P

this will highlight and select 120 records from B. hit close

4.) Right click on B, save as, "save features selected only" is ticked, and save
    as B2. this now contain just the 120 buffers

5.) Double click A, goto JOIN, and add join to B2 on rownum column added in 1 (seems to do left join still 200 results)

6.) Open attributes for A and add filter by expression "where B2_dummy is not
    null" (this makes the join effectively an inner join). This will select the 120 rows. close attribute table

7.) Right click A, save as, "save features selected only" is ticked, and save
    as A2. this now contain just the 120 polygon records from A that were within 500 metres of a point in P

Is there a simpler / better way to do this?


texts: How to make the map auto fit the available div area in full view. Is auto-repeat possible to fill white gaps

texts:
I need to show the entire world, but my map is very small at the zoom level 0 and 1. As a result it leaves white gaps and looks ugly.
Is there any way to tell the map to auto-fill the white area with map repeation ? I am using tiled layers.

Projection : 102100


texts: Transit routing using pgrouting as core routing engine?

texts:
I'm a student and I'm working on a project. The aim is given two sets of coordinates (origin & destination) route through the map to determine the optimum path using public transportation and of course walking edges.

I read about GTFS, but it's too complex, I don't care about transit schedule, time, etc. My only aim is to provide the user with a graphical representation of the shortest path, and the cost will only be the length of edges multiplied by the fare of mode of public transportation (e.g. length in kilometres * bus fare per kilometre), edge length multiplied by the average speed of the mode of transportation (e.g. length in kilometres * bus average speed by kilometre), besides, my country hasn't yet provided any GTFS data yet.

I've decide to use PostgreSQL, PostGIS and pgRouting to do this. Here's my idea and what I've done so far.

first: download an OpenStreetMap extract taken from http://metro.teczno.com/#manila.

then: convert the XML OSM data to a routable database using osm2pgrouting tool.

I've done this using Cgywin since I'm using Windows, but encountered some error though, I think? see my results below:

> reference nd=[some digits] has no corresponding node entry...
  Nodes table created
  2create ways failed:
  Types table created
  Way_tag table created
  Relations table created
  Relation_ways table created
  Classes table created
  Adding tag types and classes to database...
  Adding relations to database...
  Adding nodes to database...
  Adding ways to database...
  Creating topology...
  NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
  CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
  PL/pgSQL function "assign_vertex_id" line 14 at EXECUTE statement
  '#########################
  size of streets: 1072190
  size of splitted ways : 2281124
finished

Don't know exactly if this kind of result should concern me but I tried to do some routing. Seems working so far, if you have any comments on my osm2pgrouting result, I would appreciate the help.

Any way after converting the XML data, load QGIS application and load the ways table using DB Manager PostGIS plugin to convert the table to a shapefile, I will then enable labelling to see the edge id of each edge. then I will start tracking which edges are used for public transportation, (given that I already have a physical map which contain routes)

then: I will create a new table or multiple tables that will establish relationship between the "way" table and my public transportation data.
I haven't decided their structure yet, but that's the idea.

Now using pgRouting, I'm going to execute a shortest path query which will only route on edges that are also transit edges and walk-able by foot (footway).

Now my question is:

Am I on the right track here or am I missing something or there is something wrong with my idea?

Please I would really appreciate suggestions, corrections etc. and if there are things that would help me in making light to this project of mine.


texts: Is it possible to return 2 objects from a function when running catalog_apply?

texts:
I am normalizing a set of LAS files. In the function after normalize_height, I filter out vegetation with heights below 0 using filter_poi, and returning this as input to catalog_apply. However, I would also like to save vegetation with heights below 0 as a check.
It doesn't appear possible to return 2 objects (veg_ge_0 and veg_lt_0). It is also not possible to use opt_filter outside of the function since this seems to work on the input to the function and not on the output.
Is there a more efficient way (like returning 2 objects from the function to use in a second catalog_apply) other than creating a new function combined with catalog_apply for filtering the negatives?
library(lidR)

setwd("D:/Projects2021/p2021_0101_lidR_vs_LasTools/")
dirLASfiles_Raw <- "LasFiles__original/"
dirLASfiles_Norm <- "LasFiles_normalizedWith_lidR_wTIN/"

interpolationMethod <- tin('', extrapolate = knnidw(k = 10, p = 2, rmax = 50))

ctg <- readLAScatalog(dirLASfiles_Raw,
                      progress = TRUE,
                      filter = "-drop_withheld -drop_overlap -drop_z_below -1")

groundWaterVeg_filterAndNormalization_Function = function(chunk) {
  lasChunk = readLAS(chunk)
  if (is.empty(lasChunk)) return(NULL)

  ground <- filter_poi(lasChunk, Classification == 2L | Classification == 9L)
  if (is.empty(ground))
    return(NULL)

  vegetation <- filter_poi(lasChunk, Classification == 3L | Classification == 4L | Classification == 5L)
  if (is.empty(vegetation))
    return(NULL)

  lasNorm <- normalize_height(lasChunk, algorithm = interpolationMethod, na.rm = TRUE, use_class = c(2L, 9L))

  #lasNormFltrNegVeg <- filter_poi(lasNorm, ((Classification == 3L | Classification == 4L | Classification ==5L) & (Z < 0)), buffer == 0)
  lasNormFltr <- filter_poi(lasNorm, ((Classification == 3L | Classification == 4L | Classification ==5L) & (Z >= 0 & Z <= 30)) | (Classification == 2L | Classification == 9L), buffer == 0)  #! LASFilter is depracated. use Filter_POi
  return(lasNormFltr) #,lasNormFltrNegVeg)
}

opt_chunk_buffer(ctg) <- 25
opt_output_files(ctg) <- paste(dirLASfiles_Norm,"{ORIGINALFILENAME}_normTIN_vegge0PlusGround", sep = "")
mylasCatalogueNorm <- catalog_apply(ctg, groundWaterVeg_filterAndNormalization_Function)


texts: How to Create a Grid Inside a Polygon that is Oriented with the Longest Side?

texts:
I have been using a process called "Create Polygon Graticule" to create a grid within a polygon with assigned row and column attributes. It is a SAGA tool which I have been running in QGIS. It generates a grid based on the extent of the input polygon with the cell running N to S and E to W. It would be better for my analysis if the direction of the grid was parallel to the longest edge of the polygon. For reference I have attached images of the current automated output and an output that I manually altered to align the grid with the polygon. The grids originally covered the entire extent of the polygon but have been clipped to show only the cell parts that are inside.

Above image is what the tool out puts by default.

This is an image of a grid that was edited manually to fit the necessary parameters.

My question is, is there a way to produce grids that are aligned such as the one in the second image in an automated fashion? I will need to generate grids for a few thousand polygons.


texts: Stratified sampling with a point and polygon layer

texts:
I need to perform a stratified sample in the following context.

I have two layers; the first is a point layer (which consists of dwelling centroids.)
The second is a vector layer of polygons with each polygon being an electoral area.
There are multiple dwelling points contained within one electoral area. The dwelling point layer has all the inherited information from the electoral area layer.

I need to select 20 points from each of the five electoral areas.

How do I best do this in QGIS? I have tried all of the vector tools. - Am I missing something?


texts: Does mapCanvas().refresh() not work in QGIS 2.6?

texts:
Before I used the function

qgis.utils.iface.mapCanvas().refresh()

to reload the map canvas after for example a layer color was changed from a plugin.
But this is not working with QGIS 2.6 for me.
Do I have to use another function to refresh the map canvas or is this a bug?


texts: GEE Linear Regression Error - ReduceBands reducer must have one input:

texts:
I am having a problem with the Linear Regression function (code attached) when I try to regress latitude or longitude against B4 reflectance for a S2 image. I am getting the following error:

ReduceBands reducer must have one input: <Reducer>.

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


texts: How to modify raster pixel inside a polygon?

texts:
I would like to set to novalue all the pixels of my raster inside a polygon.

I'm trying the gdalwarp command line with -cutline option.
But this command extract the values of my raster inside the polygone (shapefile).

Raster and vector dataset's are both in the same projection.
What coud be the next steps ?
If this is the right track, is it possible to script the "raster calculator" plugins ?
Thank you in advance


texts: Programmable buffers using GDAL and python

texts: Is there a way to programmatically generate buffer layers for vector layers using GDAL and python?


texts: How to combine a vector and a raster file?

texts: I want to combine a raster and a vector file. After combining, can I query on the resulting map to obtain a specific area? How?


texts:
Is there a way to clean up shapefiles?

For example, some medians are filled in on a road shapefile I'm using.

How can I fix this kind of thing?

I'm using QGIS and working with a shapefile I downloaded from the county. In some placed the median is not filled in. I had been working with CAD in which case I can just modify the polylines.


texts: Check for NAs in `RasterLayer` and plot them in different colour

texts: I know it's an easy question. But how can I find out whats the distribution of my NA-Values to non NA-Values? And how could I plot my raster normally and make the NA values appear in a different colour than just plain white?


texts: I have locations of the points which produce sound. I want to create a sound map based of those points, topography, and the distance at which a sound can travel. On this picture you can see the red dot indicating the source of the sound; yellow dots indicate location which would be impacted by the the sound coming from the source; purple points indicate location where you won't hear the sound coming from the source location, because there is a major ridge blocking it. Is there a way to do this in ArcGIS?


texts: I'm trying to imitate the OpenTopoMap with freely available DEM-Data. Which parameters do I have to change in my hillshaded DEM to get a similar look?


texts: QGIS - WGS84 to BC Albers reprojection no curved lines, any idea why? How to snap a line to a point?

texts:
I am trying to learn a bit about QGIS.  So I thought I would drop in some layers of my Province and then recreate some points representing the BC Albers Central Meridian, and the two standard parallels.

Why when I reproject my straight lines from WGS84 to BC Albers do they not curve like they should?  I drew them in WGS84 as the image illustrates.

I hope I have outlined my problem clearly, if not please ask for clarification.

I would appreciate any input on this, thanks!

I could not post my picture illustrating this problem but hopefully it is accessible on my dropshots account: http://www.dropshots.com/stwomey%23date/2012-05-19/11:13:44

Well it looks like the link doesn't go directly to the picture...  darn it!


texts: Given latitude/longitude, how to determine whether point is in Japan or China?

texts:
For OpenStreetMap names transliteration to Latin characters, I need to know whether a point (latitude, longitude) should be transliterated to Japanese or Chinese.

What formula f(latitude, longitude) would tell whether the point is in Japan or China?

It does not need to be perfect, especially since there are islands whose status is disputed.
For transliteration purposes Taiwan is better included in the Chinese-speaking zone.

The goal is actually for the Android app OsmAnd to decide on-the-fly for thousands of labels which language to show, so it must be extremely fast, and must be able to work offline.


texts:
Hope someone can help me with this. I'm having a little trouble finding a data set which has polygons of London districts as listed here:

http://en.wikipedia.org/wiki/London_districts

I've checked Ordnance Survey and none of their boundary data sets corresponds to the above list. I have Borough polygons & ward polygons but need the above as (based on my knowledge of a couple of boroughs) these equate to sub-borough areas people would commonly refer to in everyday life.

Can anyone help? Thanks


texts: I have a dataframe of long/lat tuples and am trying to convert it into a Linestring that I can use in a Geopandas dataframe. The current format is [(-105.279, 40.060), (-105.279, 40.058), (-105.282, 40.058), (-105.282, 40.058), (-105.282, 40.03)] etc. Is there a library or technique to do this?


texts: In QGIS, the image in my print manager cannot be saved and only save one layer

texts: When I try to save an image I made from my print manager in QGIS to my computer, it only seems to save part of it. Eg one layer is saved but another layer from the image isnt. How can this be solved??


texts: Adding all feature classes of same name to mxd in ArcPy?

texts:
I have a bunch of folders which all have (almost) identical folder structures and file geodatabases of identical structures which in this folder structure. I need to add all of the same type of feature class to an mxd and merge them into a new feature class. There are a lot of them so I thought this seemed like the sort of thing I could do in python.

The folders are all in the same folder and there is a common path to the feature class in question, for example:

project_name/Subfolder1/subfolder2/project_name.gdb/featuredataset/featureclass.

so I'm thinking of trying to loop through the folders and extraxt the feature class from each of them using a method like this, but I'm not quite sure how arcpy interacts with, and loops through windows folders.

Since I'm looking for any feature class of this name within the root folder could I use a more brute force method that searches in the root folder and subfolders for the name and adds it to the mxd?


texts: Too many pixels in Earth Engine, and scale parameter meaning

texts:
I am trying to calculate stadistics for burn areas based on the number of pixels.

When running the code, GEE says too many pixels, I tried to solve it following the guides given by google developers but nithing worked. On the other hand I would like to know what the scale parameter means, for example is appropiate to use 30m for sentinel 2 imagery?

// create a function to derive extent of one burn severity class
// arguments are class number and class name
var areacount = function(cnr, name) {
 var singleMask =  classified.updateMask(classified.eq(cnr));  // mask a single class
 var stats = singleMask.reduceRegion({
  reducer: ee.Reducer.count(),               // count pixels in a single class
  geometry: area,
  scale: 50,
  maxPixels: 1e10,
  bestEffort: true
  });


texts: Distance between points inside LineString

texts:
I have a collection of geometries represented by a linestring. Each has at least two coordinates. Is it possible to write a SQL that instead of returning coordinates, return a distance between them.

input: LINESTRING( C1,C2,C3 .. CN ) size N

output: D1,D2, D3 .. DN-1 size N-1