texts: Automatically assigning attribute to line features

texts:
I have a network of railway switches and I need to assign the value of a switch line being left or right automatically.

What do I have:
A feature class containing a track that looks like this

I need to generate something that looks like the following.

The red arrow is always the direction you look at the intersection and that helps you determine if the lines branching out are left or right.

Can I achieve this in ArcGIS or using ArcPy and, if so, which library to use for coding?


texts: How do I insert a list of values into individual rows using Python?

texts:
I have a dataset that contains a field with a range of numbers (e.g. 1, 2, 3-6) as well as other fields. I've used python to expand the range field into discrete numbers (e.g. 1, 2, 3, 4, 5, 6) stored in a list in python. I would like to insert rows (6 in this example) with the range values. The values from the other fields would be copied for each new row, for each row in the table.

Row Example

| A | 1, 2, 3-6 | xxa | xxb |

to

| A | 1 | xxa | xxb |
| A | 2 | xxa | xxb |
| A | 3 | xxa | xxb |
| A | 4 | xxa | xxb |
| A | 5 | xxa | xxb |
| A | 6 | xxa | xxb |

I have adapted the code written by @Aaron. The code below works for the first row, but I would like it to read every row in the input table. In other words, it is not advancing throughout the table. I think I need to add another for loop at the beginning of insert_rows().

# import libraries

import arcpy
import os

arcpy.env.overwriteOutput = True

tbl = arcpy.GetParameterAsText(0)

def extractit(str):
    """
    Creates a list (e.g. [1,2,3,4,5,6] from string such as '1, 2, 3-6'
    """
    parts = str.replace(" ", "").split(",")

    numbers = []
    for part in parts:
        if '-' in part:
            p = part.split("-")
            numbers.extend(list(range(int(p[0]), int(p[1]) + 1)))
        else:
            numbers.append(int(part))

    return numbers

def create_table(tbl):
    """
        Creates a duplicate (but empty) table
        because we do not want to alter the original
        """
    arcpy.CreateTable_management(os.path.dirname(tbl), os.path.basename(tbl) + "_altered_v2", template=tbl)

    arcpy.DeleteRows_management(tbl + "_altered_v2")

def insert_rows(tbl):
        """
        Insert rows into new table
        """
        # Get row value in original table
        row = [row for row in arcpy.da.SearchCursor(tbl, ["MERIDIAN", "TOWNSHIP", "RANGE", "TOWNSHIP_D", "RANGE_D",
                                                          "SUBDIVIS", "NOTES_SU",
                                                          "BLMSERIA", "CATEGORY", "ACREAGE", "EXCLUSIO", "VILL_BLK",
                                                          "CONV_DAT", "IC_PAT_N", "ENCUMBRA"])][0]

        numbers = extractit(row[5])  # Extract number field

        cursor = arcpy.da.InsertCursor(tbl + "_altered_v2", ["MERIDIAN", "TOWNSHIP", "RANGE", "TOWNSHIP_D", "RANGE_D",
                                                             "SUBDIVIS", "NOTES_SU",
                                                             "BLMSERIA", "CATEGORY", "ACREAGE", "EXCLUSIO", "VILL_BLK",
                                                             "CONV_DAT", "IC_PAT_N", "ENCUMBRA"])

        for number in numbers:
            row_value = (row[0], row[1], row[2], row[3], row[4], str(number), row[6], row[7], row[8], row[9], row[10],
                         row[11], row[12], row[13], row[14])
            cursor.insertRow(row_value)

        # Delete cursor object
        del cursor

if __name__ == '__main__':
    create_table(tbl)
    insert_rows(tbl)


texts:
After unsuccessfully trying to use processing.run('qgis:idwinterpolation', ...) using a layer object in memory as the source (seems to like full file paths to physically stored datasets), I tried to use QgsIDWInterpolator instead but it crashes QGIS and creates an empty TIF file.

I'm using data from one dataset but iterating over the features to get a unique list of air quality categories and create one raster for each category. The code starts after I've already generated a temporary layer airQLayer specific to a category. I want a tif interpolation per category.  I'm not too familiar with the settings though, perhaps doing something wrong.

layer = airQLayer
layer_data = QgsInterpolator.LayerData()
layer_data.vectorLayer = layer
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute = 11  # the 12th field in the layer contains the data to interpolate
layer_data.mInputType = 1

idw_interpolator = QgsIDWInterpolator([layer_data])

export_path = rasOut # e.g. r"C:\Temp\AnalysisResults\BiomassNOx_Emissions_tpa.tif"

rect = layer.extent()
res = 10
ncols = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)
model_feedback.pushInfo(str(rect))
model_feedback.pushInfo("{} - {} - {}".format(res,ncols,nrows))
output = QgsGridFileWriter(idw_interpolator,export_path,rect,ncols,nrows)
output.writeFile()

The syntax for QgsGridFileWriter seems to have fewer parameters compared to previous examples in other threads. I thought I'd sorted it that part when it started creating a tif but no, it crashes.  Any ideas why? I'm on  QGIS 3.10.2 and have tried with and without GRASS.

Error:

QgsInterpolator::cacheBaseData qgsinterpolator.cpp:42
QgsIDWInterpolator::interpolatePoint qgsidwinterpolator.cpp:37
sipQgsIDWInterpolator::interpolatePoint sip_analysispart1.cpp:9900
QgsGridFileWriter::writeFile qgsgridfilewriter.cpp:65
meth_QgsGridFileWriter_writeFile sip_analysispart2.cpp:54
PyMethodDef_RawFastCallKeywords :
PyMethodDef_RawFastCallKeywords :
PyEval_EvalFrameDefault :
PyEval_EvalCodeWithName :
PyFunction_FastCallDict :
PyMethodDef_RawFastCallDict :
PyObject_Call :
PyInit_sip :
sipVH__core_331 sip_corepart0.cpp:17375
sipQgsProcessingAlgorithm::processAlgorithm sip_corepart2.cpp:183909
QgsProcessingAlgorithm::runPrepared qgsprocessingalgorithm.cpp:502
QgsProcessingAlgRunnerTask::run qgsprocessingalgrunnertask.cpp:66
sipQgsProcessingAlgRunnerTask::run sip_corepart1.cpp:91963
QgsTask::start qgstaskmanager.cpp:77
QThreadPoolPrivate::reset :
QThread::start :
BaseThreadInitThunk :
RtlUserThreadStart :


texts: Adding median values of several rasters (geotiff) to each fishnet polygon (200x200metres fishnet) using ArcGIS ModelBuilder?

texts:
I have a 200x200metres fishnet of my study area, with a few attributes in it, like land use, percentage of each landuse class.

Now I need to add the median value of more than 10 rasters to the fishnet, and I have to repeat this process 30 times (30 fishnets with respective +10 rasters each).

In ArcMap I know the Zonal Statistics to Table is the appropriate tool but I have not been able to produce a ModelBuilder iteration that works in giving me all the raster values within columns with their names.

I don't know how to code in Python, but if there is any tutorial on how to do this ''for dummies'' I am willing to try.


texts: How to make a virtual layer from ogr geometry?

texts:
I know how to do this:

nodes = driver.Open(node_nav,0)
node_lyr = nodes.GetLayer()
node_lyr.SetSpatialFilter(None)

But how do can use the layer methods on a geometry I created from scratch?

multip=ogr.Geometry(ogr.wkbMultiPoint)
for n in nodes:
    p=ogr.Geometry(ogr.wkbPoint)
    x,y = n.node
    x = float(x)/dim
    y = float(y)/dim

    p.AddPoint(x,y)
    multip.AddGeometry(p)

This is the traceback for this attempt:

File "/usr/lib/python2.7/dist-packages/osgeo/ogr.py", line 2927, in <lambda>
 __getattr__ = lambda self, name: _swig_getattr(self, Geometry, name)
File "/usr/lib/python2.7/dist-packages/osgeo/ogr.py", line 55, in _swig_getattr
raise AttributeError(name)
AttributeError: SetSpatialFilter

I don't want to save the geometry in a shape file (like it is described  here for example).


texts:
I collected metadata about POIs which I want to visualize. The original data however makes use of a projected CRS ("LCC projection") which I want to transform to lat/lon (WGS84):
"crs": "LCCprojection",
    "crsArg": {
        "lat_1": 46.12,
        "lat_2": 46.12,
        "center_lon": 14.815,
        "center_lat": 90,
        "x_0": -95410,
        "y_0": 6081340,
        "width": 5.67,
        "height": 4.3
    },

Unfortunately I don't have more information about the snippet above, i.e. the underlying datum etc.
I tried to transform the coordinates, but the resulting coordinates are way off:
import pyproj

lat_1=46.12
lat_2=46.12
center_lon=14.815
center_lat=90
x_0=-95410
y_0=6081340

projstr = f"+proj=lcc +lat_1={lat_1} +lat_2={lat_2} +lat_0={center_lat} +lon_0={center_lon}  +datum=WGS84 +no_defs"
print(projstr)
# +proj=lcc +lat_1=46.12 +lat_2=46.12 +lat_0=90 +lon_0=14.815  +datum=WGS84 +no_defs

source = pyproj.Proj(projstr)
print(source(x_0,y_0,inverse=True))
# (126.3370956536579, 46.67898640154669)

Since the POI is located in Maribor (Slovenia) I expect the lat/long values to be roughly around 45 and 15 degrees.
What I'm doing wrong here?


texts:
I have handwritten coordinates in the following format from a eTrex 10 device. Here are three of them.

E00197638
N00135887
E00197627
N00135939
E00197555
N00135964

They should be located in Sri Lanka, so I'm guessing the device was set up to show the UTM UPS Coordinates, he didn't write down the zone, so I'm further guessing the zone should be 44N. I have tried importing them just like that into QGIS by wriing them into a text-file (ID,x,y [Enter] 1,00197638,00135887) and using the WGS 84 datum, and UTM 44N (EPSG:32444) datum. However I arrive at random locations out of the earth or in the Indian Ocean. How can I transform these coordinates, and more importantly how can I identify the datum of those coordinates for the future.


texts:
I was trying to upload the CSV file into QGIS and somehow it doesn't read in the proper coordinate system EPSG:4326, WGS 84

Do you know what it might be the problem?? In a random exercise I tried EPSG:3857, WGS 84 pseudo mercador and the dots show up however in a completely area of the world and the point spread is rotated 90 degrees.


texts: Large scale editing to fix national park boundries in OpenStreetMap?

texts:
I recently noticed that OSM does not include the most recently updated versions of the national parks in utah since our dear leader shrank them significantly. I've taken it upon myself to fix this, and I've made up a shapefile from the state governments website. I'm fairly new to OSM editing, but from what I can tell, there isn't really a good tool for really low zoom edits - fixing boundaries this large seems hard.

What is an approach for this?


texts: How can I give random color to each line in Geoserver layer that is served from MSSQL geometry table?

texts:
I have a table of routes and I want to give each route random color to distinguish them on map (Using OpenLayers).

What will be the easiest way to do so ?


texts: Run SQL queries against shapefile in QGIS?

texts:
Is there any way I can import a shapefile into QGIS, then run spatial SQL queries against it?

I can obviously use ogr2ogr -sql to do this, but I have to import the output into QGIS each time to look at it, which is a pain.

It would be great if I could just import the shapefile into QGIS, then run SQL against it and see the results instantly.


texts: QGIS composer changing legend position when exporting

texts:
I'm observing the same issue seen in QGIS map composer changing legend position when exporting atlas as images It is happening when the checkbox "Main properties > Resize to fit contents" is checked.

Also I can observe that the Y position on page n is what should have been set as Y position calculated for page n-1.

For example is the legend was tall on page 2, and short on page 3 then it is floating high on page 3 instead of being anchored at the bottom of the page where I want it.

The upvoted answer in QGIS map composer changing legend position when exporting atlas as images won't fix the issue in my case. The only way that I found is unchecking "Resize to fit contents", which is not satisfactory for me.

Page 1, tall legend:

Page 2, short legend:

Page 3, tall legend:

I'm using QGIS 2.18.24 on OS X El Capitan.

I have the Legend properly anchored, here is the preview of page 2:


texts: Can't get ST_MapAlgebra to return fractional values

texts:
I have a 4-band raster dataset that I have imported into a Postgres/PostGIS table using raster2pgsql. I am Using ST_MapAlgebra() with a callback function to perform some band math on the dataset. I have set the return type as double precision. However, I am finding that anytime I return a decimal value (e.g., .25 or -.25) it only returns zero. However, if I return a whole number (e.g., 25 or -15), I can receive the values back from the call. I played around with this, but can’t figure out why this is occurring.

Here is the syntax of my callback function:

CREATE OR REPLACE FUNCTION urban_veg.ndvi_callback (
    value double precision[][][],
    pos integer[][],
    VARIADIC userargs text[]
)
RETURNS double precision AS
$$
    DECLARE
      r1_val double precision;
      r4_val double precision;
      result double precision;
    BEGIN
      r1_val := value[1][1][1]::double precision;
      r4_val := value[2][1][1]::double precision;
      result :=  (r4_val-r1_val)/(r4_val+r1_val);
      RETURN result;
    END;
$$
LANGUAGE 'plpgsql' IMMUTABLE COST 1000;

The call to the function is:

ST_MapAlgebra(
    ARRAY[
        ROW(rast, 1),
        ROW(rast, 4) ]::rastbandarg[],
        'urban_veg.ndvi_callback(
        double precision[],integer[],text[]
    )'::regprocedure, '16BSI'::text)


texts: Adding multiple OpenLayers Vector Layer markers to map using AJAX, JSON, and jQuery?

texts:
I am attempting to parse a JSON file and place markers on an OpenLayers based map application.  I am using jQuery to perform the AJAX and iteration on the JSON file.

I am having a problem with the code below.  When inspecting the element, I can see in the console log that the latitude and longitude are being iterated through correctly.  However, only one marker appears on the map at the 0,0 location.

Can someone point me in the right direction to make this happen?

A question similar to this was answered before (Adding multiple OpenLayers Vector Layer markers with data parsed from JSON file?), however, the jQuery based answer was replaced with pure Javascript.

 // Setup marker layer for all transmitters from DATA.json file
var markers = new OpenLayers.Layer.Markers( "Markers" );
map.addLayer(markers);

var size = new OpenLayers.Size(21,25);
var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png', size, offset);
$.getJSON('json/DATA.json', function(data) {
    $.each(data.POI, function() {
        markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(this.longitude,this.latitude),icon));
        console.log(this.longitude, this.latitude);
        markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(this.longitude, this.latitude),icon.clone()))
    });
});

Here is a sample of the JSON

"LOCATION_A": {
          "latitude": 10.2070,
          "longitude": 25.3215,
          "altitude": 1.2,
          "status": "CLOSED",
          "serverAddress": "192.168.0.1"
 }


texts:
According to the ESRI documentation, the Rebuild Indexes tool

  is used in databases and enterprise, workgroup, and desktop
  geodatabases to rebuild existing attribute or spatial indexes. In
  enterprise, workgroup, and desktop geodatabases, indexes can also be
  rebuilt on states and state_lineage geodatabase system tables and the
  delta tables of versioned datasets.

It goes on to say that the tool rebuilds the attribute and spatial indexes of base tables, delta tables, and archive tables in versioned and archive-enabled geodatabases, as well as System Tables (if the parameter is checked).  However, this and this documentation make it seem like the tool does different things if run as gdb admin vs. data owner.

In my enterprise environment, my database admin does not own data, and as such does not have options to manage the dataset via Catalog in ArcMap.

Does the geodatabase administrator login rebuild indexes for datasets not owned by the login?

EDIT:  Similarly, another esri example uses gdb admin to run the tool on system tables as well as a list of datasets (code sample below), which is what is confusing me.  I agree with what Vince is saying in the comments, but can someone verify what the tool is capable of doing?

# Execute rebuild indexes and analyze datasets
# Note: to use the "SYSTEM" option the user must be an administrator.

workspace = "Database Connections/admin.sde"

arcpy.RebuildIndexes_management(workspace, "SYSTEM", dataList, "ALL")

arcpy.AnalyzeDatasets_management(workspace, "SYSTEM", dataList, "ANALYZE_BASE", "ANALYZE_DELTA", "ANALYZE_ARCHIVE")


texts:
As part of pilot project, we have collected basic/rudimentary Household survey in one district. We have both Polygon and point data.Attributes of the database include, No of people living in the house, males, females, ownership of the house, private-public properties, structure of the house, quality of the house, etc....

I am open to suggestions on how best to visualize/ map the data (Methods). I am familiar with ArcGIS Desktop and QGIS. Recently i have started using some web data visualization tools like ArcGIS online, cartodb, Mapbox. but i am interested in Desktop solutions as i will be printing those maps for a forum.


texts: Finding and merging duplicate points in ArcGIS Desktop?

texts:
I have two point shapefiles: Red and Blue.  Each point shapefile has a z value which represents a radius (i.e. represented as buffers in the example).  A duplicate point is defined as a point that lies within Z distance of another point, which can range from 2 - 10 m.  In this example there are ten pairs of duplicate points--seven pairs overlap perfectly, while three pairs overlap by falling within another point's Z value.

I would like to keep only one of each duplicate within Z and merge the two shapefiles together in ArcGIS Desktop.

How can I find and merge duplicate points?


texts: Exporting shapefiles from large shapefile by field name?

texts:
I have a table with thousands of polyline segments. The lines in this file are named, and there are about 100 unique line names in total.

I want to create a separate shapefile for all of the segments belonging to each of the 100 unique lines.

I know I can select by attributes and export each one individually, but is there a way to automate this process?

I will be doing this many more times in the coming months and would like a better process. I am familiar with ModelBuilder but not with Python.


texts: How I can insert feature attribute after feature drawing complete

texts:
I am new to OpenLayers. I am trying to insert wfs feature information to my map after drawing like this.
I can draw a WFS  feature like (polygon, line, point) on map and after drawing the feature  complete I can insert attribute information of feature (like polygon name, area, etc) into permanently. That is after refresh the map the feature (polygon) visible on the map and on click on the feature (polygon) the feature attribute information is display by popup.
I have done drawing feature on the map by finding the some demo on Google.
But I do not understand how to store the feature attribute. Currently  I am using Geoserver  and  shapefile on  Geoserver  for layer  and  fetching  this shapefile as layer in open layer by  WFS  service of Geoserver.

Now my question is:

How to store feature attribute by using Shapefile or Database (PostgresSQL)?
If we use shapefile as WFS layer then, it is possible to store feature attribute after drawing the feature in shapefile by fetching it from Geoserver in OpenLayers.

or

For storing the feature attribute we should use Database (PostgresSQL)?
I am totally new to OpenLayers. I don’t know which way is better. I want that I am able to draw feature on  OpenLayers map  and I can insert attribute information of feature by popup and save this feature  on the map.

My website Url is: http://mindcrewgis.com/


texts: How to bring polygons closer to a point?

texts:
I have some scattered polygons with a lot of space between them.
I want to make a thematic map but all that white space doesn't look very good, so I want to bring the polygons closer to the center of the area.
Any idea on how to do that in a automatic way? It could be a raw algorithm, postgis query or python function.
Here is an example:

Original:

Final:


texts: Converting Memory Map *.qct files into QGIS?

texts:
i have created a bunch of raster base maps in Memory Map software (ie georeferenced tiff's) that are Memory map specific format (*.qct).

Is there a way to import these or convert these into QGIS at all?


texts: "Muting" first point in time animation of ArcMap

texts:
I am trying to produce a time-lapse map of points. The problem is, after enabling time, the earliest point in the layer appears before beginning the animation.  In other words, its just on the map by itself.

How can I get the first point to appear only after the animation starts, ie after I click play on the time slider?

I am using ArcMap 10.


texts:
I have a CSV with 500 ship routes.
Every row is a different trip and the columns have the values of ship_name, lat_1, lng_1, lat_2, lng_2,lat_3, lng_3... lat_20, lng_20

How can I insert this information in QGIS and create polylines between this coordinates?


texts:
I don't like the way the standard OSM maps look (too much details), and would like to check if some OSM-based services are better than Google Maps.

"List of OSM based Services" doesn't list them by popularity.

What are the main sites that use OSM maps for general purposes (ie. alternative to GoogleMaps)?

Thank you.


texts: How to maintain right-angles when buffering?

texts:
When buffering some polygons I am trying to avoid the curved buffered area like this:

I'm trying to do thins within an FME workbench but if someone knows how its done in other software that would still be helpful.

The polygons are not of a uniform shape or size so scaling them will not work. I want a uniform width of buffer. Perhaps the answer is not to change the buffer operation but to find a way of getting FME to find all curves that would be right angles and convert them into right angles. Could something like simplify do this?


texts: Removing "the_geom" field in the .csv output of Geoserver

texts:
Good day sires.

My team is currently working to develop a web-based GIS for a government client and we're using Geoserver 2.1.1 as our server application.

We're trying to play around with Geoserver's GetFeature function to generate a .csv file of the attribute fields of features saved in a PostGIS database. The function below does exactly that:

function generateFeatureTable(feature){
var featureTableURL = "http://localhost:8080/geoserver/prism/wfs?REQUEST=GetFeature&SERVICE=WFS&VERSION=1.1.0&TYPENAME=prism:" + feature + "&OUTPUTFORMAT=csv&SRSNAME=epsg:4326";
    window.open(featureTableURL);
}

The problem is that the .csv file is also including "the_geom" as one of its fields. If we make use of the PROPERTYNAME argument to specify the fields we like to include in the table, we would have to specify each and every field name excluding that of "the_geom".

Any ideas on how to specify the ones you would NOT want to include in the output file?

We will highly appreciate the insights you can provide us.


texts:
Can anyone think of a simple way to create a polygon to represent building shadows?

I have a shapefile with height information for the buildings and I want to create a simple plan looking like the image below:

I don't want it to be accurate based on the terrain etc or really specific of the sun, but just a way to demonstrate the height of the building with a shadow that is connected to the building.


texts:
I have a map viewer that was developed using the ArcGIS Runtime Java SDK. I am trying to make a print button and it works with one issue. I can't seem to figure out how to make it fit on one page.

    private void printMap() {

    Paper paper = new Paper();

    // ...set up the paper options...
    PrinterJob printJob1 = PrinterJob.getPrinterJob();
    PageFormat pageFormat = new PageFormat();
    HashPrintRequestAttributeSet set = new HashPrintRequestAttributeSet();

      if (printJob1.printDialog(set)) {
            try {printJob1.print(set);}
            catch (PrinterException exc) {
               // System.out.println(exc);
             }
         }

     //set up page format [DEFAULT] options and pass the paper settings through
    pageFormat.setPaper(paper);

    // set page format for the map to print
    map.setPageFormat(pageFormat);
    printJob1.setPageable(map);

    try {
       // send job to printer!
       printJob1.print();

       // new PrintMe().print(YourJTextComponent);

    } catch (PrinterException e) {
       // handle the exception
    }

}

I'm using Eclipse Kepler SR1, resolution is 1920x1080. Any ideas or suggestions?

This is how it prints


texts: QGIS empty date field shows value in custom form

texts: I'm trying to define custom form. Problem is with date field- form shows value even if it's empty. Is there a way to set it to show only the existing values?


texts: Creating new polygons from multiple points with centroid coordinates of polygons using ArcGIS Desktop?

texts:
How to Create from an Excel file or shapefile with multiple points with the coordinates of the centroids polygons, now I would like to create an amount x of polygons with the information of the points, it is possible to create an amount x of records for a polygon.

I have a layer of polygons of urban areas and I want to add the information of the points in a new layer of these polygons, or to replicate the polygon by each centroid point as seen in the image.


texts:
I'm working with GeoPackage files built with QGIS3.
I tried to execute some spatial queries but it was really slow, so I'd like to speed it up with indexes.

After some readings I tried

create virtual table SpatialIndex using VirtualSpatialIndex();

select * from geosirene where ROWID in (
    select ROWID from SpatialIndex where f_table_name = 'geosirene'
    and f_geometry_column = 'geom'
    and search_frame = BuildMbr(3807953,2282592,3841940,2314762,3035)
)

I tested several versions but the 'select rowid from SpatialIndex' always returns an empty set. Every table has its rtree_* . What's wrong ?

Also, I need to use geometries from other table as search_frame , is it possible ?

select b.* from buildings as b, poi_merged as p
where intersects(b.geom,p.geom) and b.ROWID in (
    select ROWID from SpatialIndex where f_table_name = 'buildings'
    and f_geometry_column = 'geom' and search_frame = p.geom)
)


texts:
I am converting osm xml to GeoJSON using ogr2ogr and I want to use the osm_id, nested under properties, as a top level id field.

Running

ogr2ogr -f GeoJSON output.geojson input.osm points

produces something like the following, where the osm_id is nested within the properties:

{
  "type": "FeatureCollection",
  "name": "points",
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
    }
  },
  "features": [
{
      "type": "Feature",
      "properties": {
        "osm_id": "32927591",
        "other_tags": "\"noref\"=>\"yes\"",
        "highway": "motorway_junction"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -122.4254065,
          37.7313829
        ]
      }
    }
  ]
}

But I would like to generate a the Geojson using the osm_id as a top level id attribute:

  "features": [
{
      "type": "Feature",
      "id": "32927591",
      "properties": {
        "osm_id": "32927591",
        "other_tags": "\"noref\"=>\"yes\"",
        "highway": "motorway_junction"
      },
      "geometry": {...

Is there a way to do this with ogr2ogr? Is there a different tool you'd recommend?


texts:
I'm trying to mask cloud cover of multi-spectral Landsat 5 image containing bands 1-7 that have been preprocessed using dark object subtraction and Landsat calibration.

I am using the fmask tool in ENVI5.5 and keeping all default parameters; Kernal Size 7 and Cloud Probability Threshold 22.5 as shown in the first PrintScreen below.

However when I try to run I get Error: Attempt to subscript list with index 0 is out of range (Second printscreen).

How do I overcome this error to produce a cloud mask?


texts: What does the prefix GP or DE on parameter data types in a Python toolbox mean?

texts:
What does GP and DE mean as seen on parameter data types.

Do they have an affect on what data is shown in the browse dialog box?


texts: Clear ALL selected features with ArcPy/Python - ArcGIS 10.6

texts:
I would like to know the Python code to clear ALL selected features. Just like the "Clear Selected Features" button works. I found this Q&A from 8 years ago:
How to Clear Selected Features in Arcmap using Python
This is the code given:
arcpy.SelectLayerByAttribute_management(Table, "CLEAR_SELECTION")
However, in the given code, you have to stipulate a table. I don't want to write code for every table/feature class. I would prefer code to clear ALL selected features. Is this possible?
Simply removing the table parameter does not work.
Edit:
Found this, which seems to be on the right track: https://community.esri.com/thread/108712
However, when I run the below code, it executes, but nothing happens:
>>> import arcpy
... mxd = arcpy.mapping.MapDocument(r"FILEPATHHERE")
... df = arcpy.mapping.ListDataFrames(mxd, "Main")
... for lyr in arcpy.mapping.ListLayers(mxd,"",df):
...     print (lyr.name)+" ...is going to be cleared now. Running 'arcpy.SelectLayerByAttribute_management'"
...     arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")
...     print (lyr.name)+" ...should be cleared now."
...
Buildings ...is going to be cleared now. Running 'arcpy.SelectLayerByAttribute_management'
Buildings ...should be cleared now.
Site ...is going to be cleared now. Running 'arcpy.SelectLayerByAttribute_management'
Site ...should be cleared now.
>>>

I ran the same thing with user2856’s suggestion:
>>> import arcpy
... mxd = arcpy.mapping.MapDocument(r"FILEPATHHERE")
... df = arcpy.mapping.ListDataFrames(mxd, "Main")
... for lyr in arcpy.mapping.ListLayers(mxd,"",df):
...     print (lyr.name)+" ...is going to be cleared now. Running...'"
...     lyr.setSelectionSet("NEW",set())
...     print (lyr.name)+" ...should be cleared now."
...
Buildings ...is going to be cleared now. Running...'
Buildings ...should be cleared now.
Site ...is going to be cleared now. Running...'
Site ...should be cleared now.
>>>

Of course, I do not actually have “FILEPATHHERE” in my code; the correct, actual file path is there (just to be clear).
All selected features remain selected when executing the above code. If I run “arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")” separately and individually on each layer, it works correctly and the selected features are deselected. What’s the problem?


texts: OpenLayers Load KML Layer in the background

texts:
I would like to show a polygon in OpenLayers. The polygon was defined in 3 KML files with difference level of detail, ie:

Low Resolution KML (with 100 points)
Medium Resolution KML (with 1000 points)
High Resolution KML (with 10000 points)

First, I created 3 vector layers for each KML file with difference max and min resolution, ie:

minRosolution: 0.0001716, maxResolution: 0.0004000
minRosolution: 0.0000429, maxResolution: 0.0001716
minRosolution: 0.0000000, maxResolution: 0.0000429

Then added those vector layers into openlayer list.
It works but OpenLayers is not usable until the whole KML is being loaded. Hence, I would like to improve it by making the loading process for Medium and High resolution KML in the background.
Is it possible?


texts:
How to set independent path in QgsOverlayAnalyzer >> intersection statement using Python?

QgsOverlayAnalyzer.intersection(QgsVectorLayer, QgsVectorLayer, QString, bool onlySelectedFeatures=False, QProgressDialog p=None)

In Qstring statement one has to define the output path and output shapefile name of the intersection. But if I use intersection inside a plugin, is there a way to set this path independent of the user? What I mean is that the path should be something like C:\users..username..etc but how can this be independent if the username? Should my plugin somehow read a path in the users computer and use this as the ouptut path? Is this possible?

My plugin uses as inputs vector layers ones that are visible in the table of contents using self.iface.legendInterface().

And another question, after the ouptut file is produced how should I open it the current map (make it visible in the table of contents)?


texts: GEE Tile error: Expected a homogeneous image collection, but an image with an incompatible band was encountered when calculating VCI from Landsat

texts:
I'm trying to calculate the VCI (NDVI - minNDVI)/(maxNDVI-minNDVI) for a collection of landsat 5 merged with a collection of landsat 8 because I will need to do this for different periods covering from 1990 to 2018. For testing, I'm using a period from 1998 to 2000 for Switzerland. If I run the following code, I run into the error:

VCI image: Tile error: Expected a homogeneous image collection, but an image with an incompatible band was encountered. Mismatched type for band 'NDVI':
Expected type: Float<-0.29947625586281806, 1.2251065502612928>.
Actual type: Float<-0.4166935481920486, 1.2483430624070735>.
Image ID: null
This band might require an explicit cast.

Is there a problem with how I calculate the NDVI? How could I solve this error?
var bbox = ee.Geometry.Polygon([
  [5.9559,45.818],
  [5.9559,47.8084],
  [10.4921,47.8084],
  [10.4921,45.818],
  [5.9559,45.818]]);

Map.centerObject(bbox);
Map.addLayer(bbox, {color: 'FF0000'}, 'geodesic polygon');

// Assign a common name to the sensor-specific bands.
var LC8_BANDS = ['B2',   'B3',    'B4',  'B5',  'B6',    'B7',    'B10']; //Landsat 8
var LC5_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7',    'B6']; //Llandsat 5
var STD_NAMES = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'temp'];

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterDate('1998-1-1', '2000-12-31')
                                                      .filterBounds(bbox);// Landsat 8
//print(l8, 'Landsat 8')

var l5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA').filterDate('1998-1-1', '2000-12-31')
                                                      .filterBounds(bbox); //Landsat 5
//print(l5, 'Landsat 5')

var withCloudiness = l8.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: bbox,
    bestEffort: true,
    scale: 30,
  });

  return image.set(cloudiness);
});

var filteredCollectionl8 = withCloudiness.filter(ee.Filter.lt('cloud', 20)).select(LC8_BANDS, STD_NAMES);

var withCloudiness = l5.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: bbox,
    bestEffort: true,
    scale: 30,
  });

  return image.set(cloudiness);
});

var filteredCollectionl5 = withCloudiness.filter(ee.Filter.lt('cloud', 20)).select(LC5_BANDS, STD_NAMES);

var lall = ee.ImageCollection(filteredCollectionl5.merge(filteredCollectionl8));
//print(lall, 'Merged')

var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['nir', 'red']).rename('NDVI');
  return image.addBands(ndvi);
};

var ndviColl = lall.map(addNDVI);

var vci = ndviColl.map(function(img){
 var id = img.id();
 var min =  img.reduceRegion(ee.Reducer.min(), bbox,300).get('NDVI');
 var max = img.reduceRegion(ee.Reducer.max(), bbox,300).get('NDVI');
 return img.expression(
   "(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   }).copyProperties(img,['system:time_start','system:time_end']);
});

print(vci, 'VCI')

// Display the result.
//Map.centerObject(vci, 9);
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
var vciParams = {min: 0, max:1}
Map.addLayer(vci.mean(), vciParams, 'VCI image');


texts: Is there a way to see the type of tags/features pulled in an OSMnx multigraph?

texts:
I'm new to OSMnx / Overpass queries.
I'm trying to pull data on rail networks in global cities. I'm using the script in Python below to pull rail data in London, for example.

There are a number of sub-categories within railway infrastructure (eg,  rail, light rail, tram, etc). I wondered if there is any way to see the kind of tags/sub categories used within the data I have pulled.

import osmnx as ox
G = ox.graph_from_point((y,x),
                        distance = 20000,
                        distance_type = 'bbox',
                        infrastructure = 'way["railway"~"rail"]',
                        network_type = 'none',
                        )

ox.plot_graph(G)

I know I could also individually go through each subcategory using 'way["railway"~"subway"]' for example, but I wondered if there was way of somehow query the actual multigraph to figure this out?


texts: How to call a method when self.iface.actionSelectRectangle().trigger() is released?

texts:
In my QGIS python plugin I use self.iface.actionSelectRectangle().trigger() to select features. I also have a button which changes the attributes of the selected features.

Is there a way instead of using a button to change the attributes, to call the method automatically after a selection was made using self.iface.actionSelectRectangle().trigger().


texts: How to find the sum of many raster files in arcmap?

texts: I have a list of 576 raster files for which I need to find the SUM. Obviously using raster calculator one by one would be tedious. Any easy script in Arcpy to do this?


texts: ArcMap calculate field based on conditional statements for other fields

texts:
I have a large dataset of species occurrence records. One field, (call it "Status_Full") contains a string of abbreviations of each designated statutory status for each record. Another field (call it "Status_Code") is a numeric field with values 1-4 donating the conservation importance of each record based on which abbreviations are in the "Status_Full" field, where 1 = "Priority", 2 = "Concern", 3 = "Locally Important" and 4 = everything else.

The "Status_Code" field has a value of 3 for any record that has an abbreviation in "Status_Full" that denotes local importance, but does not have an abbreviation that denotes higher conservation importance. At the moment a record is given a code of 3 if the "Status_Full" field contains an abbreviation that denotes local importance anywhere, e.g. if the string contains LI[VC54] ("Locally important in Vice-County 54), the record will show as locally important on the map even if it is not actually located in VC54, and hence not actually locally important. These fields were all calculated in another database before being extracted and imported into ArcMap 10.5. I now want to make a new field that gives a code value of "3" only if the previous "Status_Code" field was a 3, and the corresponding abbreviation in the string in the "Status_Full" field actually matches up with the region it is located in - contained in another field, call it "Region".

Basically I want something that says:

If [Status_Code] = 3

and [Status_Full] contains any of "abbr_1_Region_X", "abbr_2_Region_X" up to "abbr_n_Region_X"

And [Region] = "Region_X"

Then [Status_Code_New] = 3. Otherwise [Status_Code_New] = 4.

If [Status_Code] != 3, then [Status_Code_New] = [Status_Code].

If I were interested in many different regions, trying to specify all the combinations would get unmanageable pretty quickly, but I'm only interested in 3 different regions.

Below is code I tried running in the Field Calculator codeblock. I had to remove all of the link breaks to make it run. When I try to run it though, I get an "Error 99999" message, Syntax error. But all of the terms used below work to reclassify fields in isolation, it's only when I combine them all like this that it breaks. I suspect it might be easier to achieve using the Python parser, but I don't know quite where to start with that.

In the code below
[Designat_2] = "Status_Full",
[Designatio] = "Status_Code",
[Status_Code_New] = "Cat"
[NAME] = "Region"

Apologies for the confusing field names, I didn't make the database and unfortunately can't change any of them.

IF [Designatio] = "3" THEN
IF InStr([Designat_2], "CER") AND
[NAME] = "SIR CEREDIGION - CEREDIGION"
    OR
IF InStr([Designat_2], "VC46") AND
[NAME] = "SIR CEREDIGION - CEREDIGION"
OR
InStr([Designat_2], "CRM") AND
[NAME] = "SIR GAERFYRDDIN - CARMARTHENSHIRE"
    OR
InStr([Designat_2], "VC44") AND
[NAME] = "SIR GAERFYRDDIN - CARMARTHENSHIRE"
    OR
InStr([Designat_2], "PEM") AND
[NAME] = "SIR BENFRO - PEMBROKESHIRE"
OR
InStr([Designat_2], "VC45") AND
[NAME] = "SIR BENFRO - PEMBROKESHIRE"
OR
InStr([Designat_2], "WWBIC")
        THEN Cat = "3"
        ELSE Cat = "4"
        END IF
ELSE
Cat =  [Designatio]

END IF


texts: Basic User Questions: Tracing a Border in QGIS

texts:
The project: I'm plotting a dilapidated cemetery, I'm using GPS waypoint data to plot the headstones in QGIS. The map I'm making is to show the orientation of the graves, and who is in them.

I installed the Google maps plugin, now I want to trace the cemetery's road and border. Simple as that. WHY IS IT SO F** HARD TO MAKE A LINE IN THIS PROGRAM!? Online tutorials are NO HELP. I don't get why I can't just make a layer and trace the google image with a simple pen, why do I have to do a trig problem to make a line!

I uploaded the current .gpx file of the headstone waypoints and they sync up nicely though they will need some post processing. I am not experienced in GIS, but this project has to be done. I am very frustrated with this program, can you help me?

-J


texts: Getting Accurate Pixel Area in Google Earth Engine?

texts:
Is there a way to fix the pixel resolution of a dataset in Earth Engine?

I am using the USGS Mangrove layer image collection, which I reduce using the max() function. This gives me an image of 1s over the global extent of mangroves. This image was made with Landsat so its native resolution is 30m.

However, when I try to get the area of each pixel, it varies with the level of zoom I am using. If I zoom out the pixel area increases and if I zoom in, the pixel area value decreases. At the equator, the value should be approximately 900 m2 and it should change with increases in latitude.

I need to get the size of each pixel as I need to multiply it by another dataset (Hansen Forest cover) in order to get the tree cover area per pixel.

A link to my code is here: https://code.earthengine.google.com/a01b4889289544bf5cd761105d9d8caa

And a code snippet is:

Map.setCenter(9.5, 0.2, 10)

var MaxMangrove = Mangroves.max()
print('Pixel Size = ', MaxMangrove.projection().nominalScale())

Map.addLayer(MaxMangrove, {min:0,max:1}, 'Max Mangrove Value')

// not sure if it wants res in m or degrees here
var MangroveProj = MaxMangrove.setDefaultProjection('EPSG:4326', null, 30)

var PixelArea = MangroveProj.multiply(ee.Image.pixelArea())

// This runs fine and displays the image I want - an image where the
// pixel value is the pixel area. This should change with latitude
// but it changes with zoom
Map.addLayer(PixelArea, {min:0,max:100}, 'Pixel Area')


texts: Select feature in GeoJSON layer by location

texts:
I have a Leaflet slippy map with a GeoJSON layer loaded. I can click this layer manually and read attribute values from it. Now my question is:

How can I read attributes of said layer without selecting it? I'd like to do that by using my current location. I've managed to check if my location is within the GeoJSON layer with d3.js, but can't figure out how to select attributes of the respective feature. Which JavaScript client side framework could be used to face this problem?

code from d3.js console.log(d3.geoContains(jsondata[map.getCenter().lng,map.getCenter().lat]));

returns boolean answer.


texts: How to convert a raster to .bgd (ESRI binary raster)

texts:
I need to use a .bgd (ESRI binary raster) because it seems to be the only raster format accepted out of the box in the DotSpatial .NET library.

My source files are mainly .rrd/.adf/.tif files. So I instead of trying to add the gdal extension to my DotSpatial code (I tried and I failed), I am trying to convert a raster to .bgd.

My guess was that ogr2ogr should be able to do it, using a command line similar to:

ogr2ogr -f "ESRI Shapefile" C:/Temp/Shps C:/Temp/test.gdb parcel citybnd

But when I look into the GDAL raster formats I don't see anything similar to .bgd that can help me figure out which command line I should use for this .bgd conversion.


texts: Creating DEM from geoTIFF or asciiGRID LiDAR data?

texts:
For a home project I want to 3D print an area from a map.

I have downloaded a geoTIFF file and an asciiGRID file from the Tasmanian governments spatial program "TheList".

The two files are here: .asc and .tif.

I have tried for hours using separate programmes, and have found out that I need to use arcGIS, which I now have.

When I try to import the files, I come up with a grayscale flat plane.

Could anyone walk me through how to create a DEM from these LiDAR files or point me to a tutorial?


texts:
I have a .tiff image that I am trying to overlay in google earth. For that I require to georeference it first.
The image projection is UTM with NAD83 as datum. I know the lat/long of the four corners of the image.

I am doing the following way:

gdalinfo to get the four coordinates.
gdal_translate to specify lat/long info
gdal2tiles.py to tile the image.

The only thing is bugging me is that I never specified UTM and NAD83, then how come it can overlay my image?
Does anybody know how to specify this projection information?


texts: I am looking for sample data set, preferably for the US. Any map scale is okay, 1:5M 1.2M etc.  I need the data to test my reader.


texts:
Reading the solution to PostGIS > Intersect two circles. Each Circle must be built from Long/Lat degrees plus radius got me thinking.

How can I graphically display a SQL query?

It doesn't need to be cartographic quality. A sketch for debugging  purposes would be sufficient.

For example, this SQL is in the accepted solution:

SELECT
ST_Touches(temp.point1, temp.point2) as geom_touches
, ST_Intersects(temp.point1, temp.point2) as geom_intersect

FROM (

SELECT * FROM
ST_Buffer(ST_Transform(ST_GeomFromText('POINT(-105.05083 39.74823)', 4326), 2877), 1500) as point1
, ST_Buffer(ST_Transform(ST_GeomFromText('POINT(-105.04428 39.74779)', 4326), 2877), 1500) as point2

 )as temp

How could I draw ST_Touches() and ST_Intersects() - the geometry,  not the boolean?


texts: Tabulate area tool's output results units?

texts:
I'm using ArcGIS DEsktop and I have tried to calculate the area of each different land cover class in a raster of 25m size in British National Grid coordinate system. To do that I used the tabulate area tool, but I'm not sure if the results are in square meters or just meters.

Can someone help me with that?


texts: Extract Geometry information from NMEA GGA sentence in SQLite database for display in QGIS via a Virtual Layer

texts:
I have a SQLite database containing a GGA sentence in one of its fields. I'd like to display this information as a POINT layer.
I've looked at OGR Virtual Format solution as suggested here but I cannot define the Geometry field like in that solution as I do not have lat, lon readily available.
I can create these 2 fields as computed fields but I can't get them to display.

If I do this:

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>C:/Temp/SQLite/test.db</SrcDataSource>
        <SrcLayer>table_test</SrcLayer>
        <Field name="GPGGA" />
        <GeometryType>wkbPoint</GeometryType>
            <LayerSRS>WGS84</LayerSRS>
    </OGRVRTLayer>
</OGRVRTDataSource>

The GPGGA field is available and I can compute the lat and lon fields (using the Field Calculator from the layer's properties' Fields section).
However, if I add the <GeometryField> definition below to my .vrt file :

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>C:/Temp/SQLite/test.db</SrcDataSource>
        <SrcLayer>table_test</SrcLayer>
        <Field name="GPGGA" />
        <GeometryType>wkbPoint</GeometryType>
            <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="lon" y="lat"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

The GPGGA field disappear from the available fields and consequently, the lat and lon fields do not compute properly.

I have looked at the GDAL help section but couldn't find the solution.


texts: Finding correlation between point location and raster value

texts:
I have what I thought was a very simple problem but I think I may have over thought it...

I have a population density raster and a point feature layer showing the locations of X (about 70 points).

All I simply want to do is a spatial correlation; do the instances of X tend to all occur in raster cells of higher pop density? no correlation? etc. purely locational. (The points have many attributes that I could begin to do multivariate regressions on but I don't want to at this stage).

I have tried autocorrelation (Morans etc. in ArcGIS), turning X into a Boolean raster and running band statistics in ArcGIS, Geographically Weighted Regression and I'm now onto Bivariate Spatial Association.

Someone tell me I've had a massive oversight?! Perhaps plotting the raster cell values against a Boolean 1 or 0 of the variable X?

This question is similar but fizzled out while this person looked like they want to test an interpolation. Don't mind answers in Arc and/or R.


texts: Creating checkboxes to select multiple inputs using ArcPy?

texts:
I am working on a script (to publish as a geoprocessing service) that will take multiple feature layers and use them as inputs to the ExportCAD_Conversion tool.

What I ultimately want is a set of checkboxes that allow the user to check which layers they want converted into CAD (from a list of about 10 features).

I saw somebody mention using something like this:

import arcpy

Feature = "C:\FeatureLocation"

    FeatureCheck = arcpy.GetParameterAsText(0)

    if str(FeatureCheck) == "C:\FeatureLocation":
            arcpy.AddMessage("Feature")
            result = arcpy.Exists(Feature)

    Feature2 = "C:\Feature2Location"

    Feature2Check = arcpy.GetParameterAsText(1)

    if str(Feature2Check) == "C:\Feature2Location":
            arcpy.AddMessage("Feature2")
            result = arcpy.Exists(Feature2)

    arcpy.ExportCAD_conversion([FeatureCheck,Feature2Check], "DWG_R2010", r"C:\dwg\test.dwg", "Ignore_Filenames_in_Tables", "Append_To_Existing_Files", "")

This only got me an error that said "Dataset true does not exist or is not supported"

Am I going down the right path?


texts: Navigation zoom in / zoom out not working in ArcGIS API for JavaScript?

texts:
is it just me or the zoom in and zoom out customized buttons [The ones at the top with the icons that read "Zoom in" and "Zoom out"]do not work?

http://developers.arcgis.com/javascript/samples/toolbar_navigation/

[ Regular zooms (the '+' and the '-') do work. ]

I'm currently developing a gis application using arcgis 10.2 javascript api and can't get those buttons to work.


texts: Truncate coordinates (LonLat) to decimal precision

texts:
For example if coordinates are 17.12345678, -82.12345678, and I want to convert them to 3 decimal digits.

Expected result: 17.123, -82.123

Is there any API in OpenLayers 2 to truncate coordinates (LonLat) to decimal precision?


texts: In the esri geodatabase, what is the difference between Simple and Attributed relationship classes?


texts: What causes this error when loading plugin

texts:
After installing pyspatialite on my Mac two plugins fail to load with this error::
This plugin is broken
invalid literal for int() with base 10: '2-spatialite'

I am fairly sure that the offending string '2-spatialite' is part of the version string associated with pyspatialite (2.6.2-spatialite.2.4.0-4)

The two plugins giving this error are Qspatialite and opengeo-explorer. What is weird is that the opengeo plugin does not use spatialite.

I assume that something in the plugin support is splitting the version string on '.' and expecting to get integers back.


texts:
If "enable on the fly" is on and the CRS changes at all during the lifetime of a project then seemingly, it refuses to measure (using the measure tool I just get 0.00000). So new project with all the layers added in the same CRS (OSGB British National Grid EPSG:27700 which we used one QGIS 2.18) it measures, but it gets it wrong. When compared with Google it's about 0.8m rather than 1m.

If the project is all saved as and only ever goes into EPSG:3857 then it seems to be right, but if you accidentally add a layer with a different CRS, it doesn't measure at all.

Am I doing something wrong?


texts: Defining parameters for Python Toolbox tool in ArcGIS Pro

texts:
I have been posting similar questions recently, but I found a template code for a Python toolbox and I would like more help.

I am trying to create a tool in ArcGIS Pro that takes parameters from users (initially, the geodatabase filename) and shows the report on the screen.

I added to the template the function "getParameter." I thought this function will open up a dialog box for a user to choose a parameter, but it does not for some reason.

The following code has no syntax error. When I add a new Python toolbox and ran this code, it does run, but it says "no parameters," as the screenshot says. What change do I need to make?

import arcpy
import os

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

    # List of tool classes associated with this toolbox
        self.tools = [Report]

class Report(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Tool"
        self.description = ""
        self.canRunInBackground = False

    def getParameter(self):
        arcpy.env.workspace = arcpy.GetParameterAsText(0)

        featureclasses = arcpy.ListFeatureClasses()
        for fc = in featureclasses:
            print (fc)

    def getParameterInfo(self):
        """Define parameter definitions"""
        params = None
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

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

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

    def execute(self, parameters, messages):
        """The source code of the tool."""
        return


texts: Order of parameters in QGIS 3.6 graphical modeller

texts:
When using graphical modeler in QGIS 3.6 I have added additional parameters to the model inputs for the buffer tool.

The model looks like this:

When running the tool the order of the parameters in the resultant window do not follow the order parameters were added. They appear to be listed alphabetically:

I can change this quite easily by exporting the model as a Python script (as per this question) and changing the order of the input parameters there:

This means my model is a now a Python script rather than a .model3 file. Is there any way to change the order of parameters through graphical modeler. Or alternatively is there a way to save a Python script as a .model3 file?


texts: How to set polgon fill type to no color in Mapbox

texts:
I'm new to Mapbox and I am trying to remove any color in my polygon fill using the Mapbox API. Note: I don't want to set opacity/transparency. I simply want there to be no color. It keeps defaulting to "#000000", even if I comment out the 'fill-color' property. Is there any way to achieve this?
 map.addLayer({
            'id': 'us-states',
            'type': 'fill',
            'source': 'us-states',
            'layout': {},
            'paint': {
                //'fill-color': '#333',
                'fill-opacity': 0.8,
                'fill-outline-color': '#000'
            }
    })


texts: Unable to submit changes in feature layer using applyEdits arcgis for android

texts:
After fetched my webmap from portal and found the layer I want add a polyline to, I'm trying to use the applyEdits method to update it in portal. Unfortunately, it doesn't work.

My code goes as follows:

    Pesticide pesticide = globals.getPesticide();
    Map <String,Object> attributes = new HashMap<String, Object>();

    attributes.put("Name",pesticide.getAgricName());
    attributes.put("FarmeName",pesticide.getSiteName());
    attributes.put("CropType",pesticide.getCropType());
    attributes.put("ParcelID ", pesticide.getSiteID());
    attributes.put("PesticideType ", pesticide.getPestType());

    pestLayer = globals.getWebMapLayer();
    graphic = new Graphic (polyline, lineSymbol,attributes,0);

   Graphic [] adds = {graphic};

    ArcGISFeatureLayer featureLayer = new ArcGISFeatureLayer(pestLayer.getUrl(), ArcGISFeatureLayer.MODE.SNAPSHOT);
   if(featureLayer.isAllowGeometryUpdates()) {
       Toast.makeText(getApplicationContext(), "all geometries are able to update", Toast.LENGTH_LONG).show();
   }
   else{
     Toast.makeText(getApplicationContext(), "all geometries are NOTTTTTTT able to update", Toast.LENGTH_LONG).show();
   }

    featureLayer.applyEdits(adds, null, null, new CallbackListener<FeatureEditResult[][]>() {

        @Override
        public void onCallback(FeatureEditResult[][] featureEditResults) {
            if (featureEditResults[0] != null && featureEditResults[0][0] != null && featureEditResults[0][0].isSuccess()) {
                Toast.makeText(getApplicationContext(), "Is Success", Toast.LENGTH_LONG).show();
            }
        }

        @Override
        public void onError(Throwable throwable) {
            Toast.makeText(getApplicationContext(), "failed", Toast.LENGTH_LONG).show();
        }
    });

Can anyone help me?


texts:
I have a fairly small script which is running a cursor on a small feature class (~2500 features) which is crashing python with a python.exe has stopped working message.

I have a cursor that loops over features, deletes any that don't meet a set criteria, looks up field value keys in a dictionary, and updates another field with the dictionary value.  This works well.  But if I add in a check to deleteRow() any records that have a Null or a value that doesn't start with the letter O.

arcpy.FeatureClassToFeatureClass_conversion(r"mydatabase.sde\mydatabase.dbo.SitePolygon", "in_memory", "xTempPolygon")

with arcpy.da.UpdateCursor(r'in_memory\xTempPolygon', ['RecID', 'LotID', 'Owner', 'Manager', 'Class']) as cursor:
    for row in cursor:
        if row[2] != 'CITY' or row[4][:4].upper() == 'ROAD' or row[3].upper() in ['WASTEWATER', 'STORMWATER', 'WATERSUPPLY']:
            cursor.deleteRow()
        elif row[0] in resultDict:
            row[1] = resultDict[row[0]]
            cursor.updateRow(row)

        #################################################
        # If I add this block, it crashes
        if row[1] is None or row[1][0].upper() != 'O':
            cursor.deleteRow()
        #################################################

The resultDict is a dictionary that comes from a SQL lookup that returns RecID and LotID.  These are looked up in the cursor and populates the LotID on the feature based on its RecID.

If I leave out the final if block, it runs without issue.  If I add it in it crashes.  I need it there because I'm ending up with some records that have a NULL or a value that doesn't start with O in the LotID field.  I need these records deleted.

I need the extra if to run after the update to allow for the update to run, as no records have any values in row[1] before the update - they are all NULL so they'd all be deleted.

What is causing this crash, and how do I avoid it?


texts:
I'm trying to import a large dataset to Postgres, of about 600MB, that sits on a remote AWS server. But I get the following error

Not enough memory to allocate requested memory (nRecords=26923892).
Probably broken SHP file

I have a single core, 2GB ram machine. My question is, should I scale up the machine or if only I increase the memory allocation I get shp2pgsql to run smoothly even for big datasets (600MB and bigger)


texts: Create a regular grid of sample points with angle and give them a value

texts:
I'm mapping able points to archaeological research.
Using QGIS, I need to create a regular grid of sample points, give them a value, specifying spacing and rotation. Better if it results inside a polygon.

I'd find some plugins and core actions in QGIS but none of them fully achieved my goal.
ArcMap has Create Fishnet that is quite close what I'm looking for but it's not possible to create fields with incremental value allowing points to have ID (A1, A2... B1, B2...). After create it, as Create Fishnet output results in line or polygon feature. I have to use Intersect and Merge to create a point grid.

I'm a beginner in programming. Maybe I can have some tips how to start focused on that.

I can describe a process:

Choose origin grid reference - Maybe a single point feature on a polygon corner or coordinate
Input an angle value - Measured in one side of the polygon or opposite corner
Space between point in x-axis
Space between point in y-axis
Number of columns x-axis
Number of lines y-axis
Create field X_AXIS: Input incremental string value by row(A,B,C...Z,AA,BB,CC...)
Create field Y_AXIS: Input incremental integer value column (1,2,3...)
Create field POINT_ID: Concatenate X_AXIS and Y_AXIS


texts: QGIS 2.4 plugin PyQT qWebview javascript interface problem

texts:
I have a QGIS plugin which contains a qtwebview and Javascript communication between the page and the plugin.

It is working perfectly with QGIS 2.2  but not with 2.4  (both QGIS versions share the same plugin directory),

The javascript side never gets called with evaluateJavaScript() and the python side never gets invoked from Javascript either.

What could explain that it is working in 2.2 but not in 2.4?

def on_loadFinished(self):

...

def init_Clicked(self,iface):
    self.ui.webViewResult.page().mainFrame().addToJavaScriptWindowObject("QGISInterface", self.moTAZJsInterface)
    self.ui.webViewResult.loadFinished.connect(self.on_loadFinished)
    self.ui.webViewResult.load( QUrl(config["TAZ_URL_INTERFACE"]))
    self.ui.webViewResult.show();


texts: Georeference vector based on same data as raster?

texts:
It sounds complicated, but this approach would be quite simple. I'm just not sure, if I can solve my problem in this easy way.

I have a raster image and a set of vector data, that are all not georeferenced and stored in a vector program. The task is to georeference all data, but keep the image as raster, and the vector data as shapefiles.

The plain raster image I have already georeferenced.

Now I am wondering, if the following approach would work:
I could take the raster image with the vector data on top and export that all together as image file. As I have the raster already georeferenced, georeferencing this raster should be easily, as it is exactly the same, only with a few more data (plain image information) on it.
Could I now take each vector dataset and georeference it based on the already georeferenced raster with the vectors (as image information)? This would result in setting the links between the real vector data (the one that should be spatially adjusted) and the same vector data (but as image information). Spatially adjusting would in this case just be to always click on the same point between the real vector and the same vector, which is on the raster though.

I am not sure if my explanation makes sense, or how I could better describe my idea. Hopefully anyone understands me.


texts: Nearest point of MultiPoint to another Point, over a LineString

texts:
I need to calculate the distance between some points on a line (LineString), using PostGIS, to find the nearest over a set of MultiPoint.

I know there is a function ST_Distance(), but this one returns a distance as the crow flies, with no possibility to take a path into account.

The problem is : "which of A or B is the nearest of S ?"

Here is a little drawing, for a better explanation :

Is there an alternative to ST_Distance() ?

I have already read Getting distance between two points over linestring using PostGIS?, but unfortunately, it doesn't give a complete reply : I have a MultiPoint feature.

The solution given wasn't succesfull, but I guess it may just need an adaptation ..?


texts: How to print log from FME when run from python script?

texts:
I created a Spatial ETL ToolBox in ArcCatalog.
When I run the tool from ArcCatalog I see log massages.
But when I run this tool from a python script, the logs do not print to the python console.

import arcpy
arcpy.ImportToolbox(r"c:\Toolbox123.tbx")
arcpy.SpatialETLTool(param1,param2,...)

How to print log from FME when running from python script?


texts: Use Function V.clean of the grass tool in a Python script

texts:
I need to integrate V.clean functionality of Grass in a script python. this functionnality is accessible in Qgis (Processing->tools).
I found that I should use Grass API but I don't know the code syntax, Please give me examples of use of the Grass API functionalities in a Python script, especially the v.clean function


texts: Link Leaflet opacity slider to uppermost layer

texts:
From the Leaflet.OpacityControls plugin, you can add an opacity slider to a tiled layer by:

var opacitySlider = new L.Control.opacitySlider();
map.addControl(opacitySlider);
opacitySlider.setOpacityLayer(interactive_layer);

I want the opacity slide-bar to link to not just one layer, but the tiled raster layer with the largest z-index (i.e. upper most layer) at any given time. How should I replace the interactive_layer var to do so?
Working with very limited JavaScript experience.


texts: Working example to group layers in LayerSwitcher of OpenLayers 2?

texts: I am looking for a nice plugin with working examples to group layers in OL 2.11/2.13.1 LayerSwitcher. I found several emails and posts explainig that GeoExt is the best way to do this, but have only access to OL.


texts: Adding only unique features from another shapefile

texts:
I am new to a company that makes use of a rather large database of point data. For whatever reason, the woman who was here before me made a huge amount of duplicate points - the attributes may be a little different but the coordinates are the same.

I've found that I can remove the points either in QGIS using the mmqgis plugin, or in ArcMap using the "delete identical" tool. However, when I use either method, I have noticed that there are some points that are removed completely. For most, it leaves one point but I have found a few that are completely gone.

Do any of you know of a method to add unique features to an attribute table?

I'm thinking that I can compare the original table to the "cleaned" one and add the missing points back. Loo


texts: How to map in Graduated style in QGIS when it does not show the option to?

texts: I am trying to map a column through Graduated style in QGIS. But the values in my column have a few repeated values of example 1000s, 1003s, 1005s etc. and the program recognizes it as categorized style. If I map it as categorized I will have over 1400 categories. What can I do to plot as a Graduated style? It is not appearing in the drop down option.


texts: Reduce the resolution of a shapefile in R (or QGIS if necessary)

texts: I'm using shapefiles of history U.S. state boundaries from NHGIS to make maps in R. These shapefiles are much higher resolution than I need when making maps at the scale of the country: over 2 million observations for one file. Is there a way that I can reduce the resolution of these shapefiles in R? If necessary, I can do it in QGIS instead.


texts: Creating chloropleth maps in geopands with 'user_defined' bins

texts:
As suggested in the geopandas mapping tools documentation, I want to create a chloropleth map using a custom defined bin sizes:

world.plot(column='gdp_per_cap', cmap='OrRd', scheme='user_defined')

Indeed, 'user_defined' is mentioned explicitely in the geopandas plotting.py code on github.

However, when I attempt this, python (v 3.5, geopandas v 0.4.0) throws an error:

ValueError: Invalid scheme. Scheme must be in the set: dict_keys(['quantiles', 'equal_interval', 'fisher_jenks'])

Indeed, when inspecting the source code in my environment, the schemes are restricted to just these four.

How can it be that my Geopandas installation is not in sync with what's being suggested in the Geopandas docs and source code? According to Conda, I really have the latest version of Geopandas installed.

ANSWER:

I've cloned the Geopandas repo from Github. Using the plotting functionality from there, the following works:

world.plot(column='gdp_per_cap', cmap='OrRd', scheme='user_defined', classification_kwds={'bins':[1000, 10000, 100000]})


texts: Problem with PostgreSQL/QGIS converting long/lat to suitable projection/SRID

texts:
I'm sure it's a simple question but this is all very new to me and after a lot of internet searching my points are still in the sea!

I have a set of UK postcode data that I have converted into latitude and longitude. I also have a UK counties shapefile.

The code that I have used to convert the latitude/longitude to a geometry point column is as follows:

ALTER TABLE colic ADD COLUMN longlat geometry(POINT,27700);
update colic set longlat = st_setsrid(st_point( longitude, latitude), 27700);

And this produces a column that has values as follows (I've only included a couple!):

  "0101000020346C000048BF7D1D3867E83FC05B2041F1334A40"
  "0101000020346C0000F085C954C1A8F7BFA1F831E6AE954A40"
  "0101000020346C000020D26F5F07CEF4BFE3361AC05B504A40"
  "0101000020346C00007D3F355EBA49F73FB6847CD0B3614A40"

The .prj from my shapefile is:

  PROJCS["British_National_Grid", GEOGCS["GCS_OSGB_1936",
  DATUM["D_OSGB_1936", SPHEROID["Airy_1830",6377563.396,299.3249646]],
  PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["False_Easting",400000.0],
  PARAMETER["False_Northing",-100000.0],
  PARAMETER["Central_Meridian",-2.0],
  PARAMETER["Scale_Factor",0.9996012717],
  PARAMETER["Latitude_Of_Origin",49.0], UNIT["Meter",1.0]]

And I set its SRID to 27700 when I upload it using the PostGIS shapefile uploader.

In PG Admin, when I ask:

select st_srid(geom) from counties limit 1;
select st_srid(longlat) from colic limit 1;

I get 27700 for both.

But... when I try any spatial query I get no response, and when I load the layers into QGIS, the point data are off the tip of Cornwall in the North Atlantic, and only represented as one dot on the scale of the counties map.

If I zoom in on them, they are 'there', but just not in East Anglia where they should be!

This is my first attempt at both SQL and mapping - I'm sure there is something really simple that I have missed but I would be so grateful if someone could help me.


texts: ArcMap Image Analysis Window - functionality is grayed out

texts:
In ArcMap's Image Analysis Window, none of the tools are available even though I have a raster selected in the window.  I do not have Editor turned on.  I can't identify why the Processing tools are not available to me.

Any ideas?


texts: Query feature class on attribute value, and generate Python list of feature classes

texts:
I want to query a feature class on attribute 'state' == 'new', and if so, add that feature class record to a list (I will process this list later).

I have:

cursor = arcpy.da.SearchCursor(points, 'state')
for row in cursor:
    if row[0] == 'new':
        new_point_list.append(row)

('points' is the feature class, 'state' is the attribute within the feature class)
but that just makes a list of tuples, not feature class records!

I ultimately want a list of records in that feature class with a status of 'new', so that I can loop through them and perform other functions.


texts:
I am developing a web application that dynamically creates tables of data on CartoDB. Once I've uploaded some data, I'd like to show a cluster of points from that data. The uploaded data has a the_geom field of datatype geometry with well-identified data (eg, ST_SetSRID(ST_Point(37.878086,-122.51026),4326))

So, given that I get a table name back from a successful create/write (eg, t62b61b0091e633029be9332b5f20bf74_6d6d454828c05e8ceea03c99cc5f5), I'd like to grab the viz.json identifier (eg, 2b13c956-e7c1-11e2-806b-5404a6a683d5, or even the full https://myCartoAccount.cartodb.com/api/v2/viz/2b13c956-e7c1-11e2-806b-5404a6a683d5/viz.json) somehow, programmatically.

I can't seem to find the way to do this, nor references to it on the documentation.

Help?


texts: Multiple conditions for new field return NULL value in QGIS when using field calculator

texts:
I have a field with building height in levels.
I want to create a new field that returns 'low-rise' if buildinglevels < 4. That returns 'mid-rise' if building level is 4-9 and that returns 'high-rise' if building level > 9.
I can do this for one condition

If( "buildingle" > 9, 'high-rise' ,pass)

But when trying the others, they overwrite the previous to NULL.
How to solve this?


texts:
I currently have a wonderful mobile version of our site up in which it uses PopupMobile to display the requested info. I also map a button on the map that allows the user to do a random function.

The problem that I am running in to is since the button is on top of the map, when a user expands the mobile popup the button is on top of it as well. Is there a way to keep the button on top of the map, but make it go underneath or not be there when the popup is expanded?

I have included a small picture illustrating my issue. The left side is when there is just a map there, and the right is when a user expanded the popup infowindow.


texts: The directions in http://docs.qgis.org/2.2/en/docs/user_manual/osm/openstreetmap.html#saving-osm-data seem to have been written before the OSM plugin was merged into QGIS proper. I can locate no OSM plugin or toolbar, and the vector menu only has options for opening a .osm file in QGIS.


texts: Exporting map with arcobjects interupted by refresh event

texts:
I have developed an ArcMap extension in C# for ArcGIS 10.0 with a Form that configures the PageLayout and exports the ActiveView as PDF. There are two different ways to do that. You could use IPrintAndExport.Export function or you could use the IActiveView.Output function. Both ways generally work fine and export a PDF document as expected as long as the user is patient enough to wait for the export to finish without clicking somewhere or moving the dialog. If the user can’t keep his finger still during the exporting, there seems to be thrown an event that forces a refresh on the ActiveView, immediately interrupts the map export and generates half-finished pdf. I just tried to prevent the user to do something by implementing an IProgressDialog2 but again with the same result. When the user moves the IProgressDialog2 map export stops again.

When I use the ESRI GUI standard function accessible through the main menu ‘file --> export map’ the whole ArcMap application is freezed during the export. So every click of the user is ignored and the map export is finished successfully. Any idea how to do this in arcobjects? Or is there a different solution?

Attached you will find my source code:

string exportFileName = @"C:\test.pdf";
double screenresolution = 96.0;
IExportPDF pdfExport = new ExportPDFClass();
IExport export = (IExport)pdfExport;
export.Resolution = 300;
export.ExportFileName = exportFileName;

ESRI.ArcGIS.esriSystem.tagRECT exportFrame = ((IActiveView)m_mxDoc.ActiveView).ExportFrame;
double right = exportFrame.right * export.Resolution / screenresolution;
double bottom = exportFrame.bottom * export.Resolution / screenresolution;
exportFrame.right = (int)(right);
exportFrame.bottom = (int)(bottom);
IEnvelope pDisplayEnv = new EnvelopeClass();
pDisplayEnv.PutCoords(exportFrame.left, exportFrame.top, exportFrame.right, exportFrame.bottom);
export.PixelBounds = pDisplayEnv;

Cursor.Current = Cursors.WaitCursor;

int hdc = export.StartExporting();
m_mxDoc.ActiveView.Output(hdc, (int) export.Resolution, exportFrame, null, null);
export.FinishExporting();
export.Cleanup();

Cursor.Current = Cursors.Default;


texts: Set order of polygons to burn into a raster with GDAL?

texts:
I have intersecting polygons with an attribute that I'd like to burn into a raster with GDAL. In this simple example, the attribute Value from three polygons in one shapefile is to be transferred to a one-band raster. However, I have hundreds of polygons and they are not organised according to attribute values or location.

If I run gdal_rasterize:

  gdal_rasterize -ot Float32 -of gtiff -init -a Value -tr 1 1 -te -150 -150 150 150 -l test  test.shp test.tif

I get a nice raster:

but my problem is that I want the highest value of any intersecting polygon to be burnt into the raster, not the last polygon.

What is the easiest way to burn the highest attribute value into a raster, only using gdal/ogr?


texts: Conversion of any projection to latlong using proj4 string

texts:
I am using a tool which uses proj4 lib to project maps.

I have a set of xy coordinates to be converted to latlongs with that tool. as I have the source projection string say

+proj=aea +datum=NAD83 +ellps=GRS80 +lon_0=-96d +lat_1=45d30 +lat_2=29d30 +lat_0=37d30 +x_0=000000 +y_0=00000

if I pass

+proj=longlat +datum=NAD83 +ellps=GRS80 +lon_0=-96d +lat_1=45d30 +lat_2=29d30 +lat_0=37d30 +x_0=000000 +y_0=00000

as the destination projection

Will it give me a LAT LONG from the x,y supplied?

I know there are ways to write code and do inverse projection to get latlong, but I wanted to bypass that and use x,y as the input to the tool I have and get latlong from it by just using the proj4 string.


texts: QGIS and the South Pole: bug or feature?

texts:
I am trying to work with some data from Antarctica in QGIS. I have defined this new projection:

+proj=laea +lat_0=-90 +lon_0=0 +datum=WGS84 +units=m +no_defs

And started adding a few vector layers. One of these contains the location of meteorologic stations in the continent. No matter what I try, QGis never portraits the South Pole station:

Is this the result of an issue with the custom CRS or is it a bug with QGis?


texts:
I saw this question severall times (e.g. What projection to use to get a game map into QGIS? and How do I specify a CRS for a fictional game map in QGIS?) but since I am new to GIS, I need further instructions.

I am trying to port the map from the game PLAYERUNKNOWN's BATTLEGROUNDS to the right format for a tileserver. My goal is to create an android app which shows the map using some SDK like Google Maps or Mapbox. There already exists an online version from another developer.

So far I used MapTiler to convert an image of the map to map tiles. However, when I use "Standard Tiles", the map is projected on a globe an therefore the image is stretched and clinched. Using the option "Raster Tiles" creates some white borders around my map and I think makes it harder to add POIs later, because I don't have a coordinate system.

Can anyone tell me how to convert this image of the map into a GIS format while still preserving the coordinate system from the game so that I can later add POIs (spawn locations for vehicles and items) on the map using the same coordinates as in the game?

For reference: The distance between two white lines on the linked map is 100m and the distance between two yellow lines is 1km.

I am an experienced android developer but this projections and coordinates stuff is very confusing for an outsider.


texts: How to transfer data from shapefiles into a PostGIS database?

texts:
I have a project in QGIS that has 15 vector layers, each of which is from an ESRI shapefile, and represents the property parcels in a county in Washington State. Each file has an attribute table with different field names (for example, one has a field called owner and another has two fields for the same data called owner_first and owner_last).

I want to move this data all into one PostGIS database with the data normalized into consistent fields, with an additional county field to identify which county the data was from. Are there any QGIS built-ins or plugins for moving data from shapefiles into PostGIS, and doing whatever conversions are necessary to re-map the shapefile attribute table fields to the database schema?

I understand that it will be slightly different for each shapefile, based on the way that particular shapefile is structured, but I'm just trying to get a handle on the basic process. As a simple example to illustrate what I'm trying to do, suppose my PostGIS database was called parcel_data and it had a table called parcels with the following fields:

id         (Integer)   # primary key
parcel_no  (String)    # parcel number (string bc can be of form 0000123433)
owner_id   (Integer)   # a foreign key into owners table
address    (String)    # street address
value      (Float)     # property value
county     (String)    # name of county

and another table for owners:

id      (Integer)      #primary key
name    (string)       #owner's full name
address (string)       #street address
zip     (string)       #zip code

and I wanted to import data into these tables from a shapefile called mason_county.shp that had the following fields in its attribute table:

parcel           # insert into the parcel_no field in DB
owner_last       # last name of owner
owner_first      # first name ... concatenate w/ last & insert into "name"
street_addr      # the street address, inserted into "address" field
propvalu         # should be mapped to "value" field in DB
zip              # zip code - insert into "zip field"

So I want to insert all of these into the correct DB tables, along with the geometry/CRS/etc information, and then use "mason" as the value for the county field in the parcels table. How would I go about doing this? I am using QGIS, and can also use tools like ogr2ogr.  Is there a way to do this from inside of QGIS? If there is no way to do it with QGIS directly or from command line with ogr2ogr, I am also open to suggestions using Python.


texts: Finding the origin in .laz file with LAStools

texts:
I obtained a .laz file, which I used lasinfo and las2txt to view the records.

 Scale Factor X Y Z:          0.00100000000000 0.00100000000000 0.00100000000000
  Offset X Y Z:                22717.000 31027.000 21.000
  Min X Y Z:                   23776.000 30560.000 14.859
  Max X Y Z:                   23807.997 30591.999 43.696
  Spatial Reference:
PROJCS["SVY21 / Singapore TM",
    GEOGCS["WGS84, WGS84",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",1.366666666666667],
    PARAMETER["central_meridian",103.8333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",28001.642],
    PARAMETER["false_northing",38744.572],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

and a sample record is 23776.000,30560.047,19.450.

Now, how do I interpret the above information, as in where is the origin in which the record referenced? When I punch in the lat, long from the above, I got the point below. If that is the origin and the X and Y is pointing eastwards and northwards respectively, then my assumption is incorrect, as I am certain that the data originated from somewhere southwest of that marker.

.


texts: How can I visualize each prior version of a OSM way?

texts:
What tools are available to visualize each historic version of an OSM way?
In this case I want to find out the history of the (incorrect) long straight sections on Eagle View Trail and Vista Grande Road:

https://www.openstreetmap.org/way/127011426

How can I visualize each version and then track the matching change sets?


texts: Force network analyst route to follow a certain pattern

texts:
I'm using ArcGis 10.4.1 and network analyst. I'm trying to force,  network analyst route to follow a certain pattern.

I've tried editing network attributes in network dataset properties and changing values in Global Turn Delay Evaluator, but with no success.
Instead what network analyst is shown below.

.

Can anyone suggest me any solution. I believe that i have to hard code the route behavior.


texts:
I am using the Esri.ArcGISRuntime 100.5.0 NuGet package with C#.

I want to add something like a border around a point symbol. So I get the geometry of the symbol (this can be a circle, square, etc.) which is a Polygon geometry.

Based on this geometry I want to create a simple red Polyline.

I tried the following, but I only get a filled rectangle:

var solidStrokeSymbolLayer = new SolidStrokeSymbolLayer(solidStrokeSymbolLayer.Width, Color.Red)
{
    IsColorLocked = true,
    IsEnabled = true
};
var polygonSymbol = new MultilayerPolygonSymbol(new List<SymbolLayer>
{
    solidStrokeSymbolLayer
});
var vectorMarkerSymbolElements = new List<VectorMarkerSymbolElement>
{
    new VectorMarkerSymbolElement(geometry, polygonSymbol)
};
var vectorLayer = new VectorMarkerSymbolLayer(vectorMarkerSymbolElements)
{
    IsEnabled = true,
    IsColorLocked = true
};

var pointSymbol = new MultilayerPointSymbol(new List<SymbolLayer>
{
    vectorLayer
})
{
    Size = 30d
};


texts: Installation of Whitebox plugin in QGIS 2.18

texts:
I want to install the plugin Whitebox (https://plugins.bruy.me/processing-whitebox.html) into QGIS 2.18 but it is not listed in the Plugin Repository. I have the zip file. Is it possible?

I have given up on QGIS 3.xx for my current tasks for the time (where I can install the plugin).


texts: Canadian province boundaries in latitude and longitude

texts: Are there any datasets available that give a list of latitude and longitude points for Canadian province boundaries?  I don't need anything too detailed as I'm displaying collated data in smallish info graphics.


texts: Counting raster values using another raster as mask in R?

texts:
I would like to count the frequency of values in one raster, using another raster (different resolution) as the zonal condition.

The premise is that the coarser resolution raster represents the area that I would like to count the frequency of the finer resolution raster;

# make 2 rasters of the same extent, different resolutions
ext <- extent(0,1000,0,1000)
r1 <- raster(nrows=1000, ncols=1000,ext)
r1[] <- sample(seq(from = 1, to = 6, by = 1), size = 1000000, replace = TRUE)
r2 <- raster(nrows=10, ncols=10,ext)
r2[] <- sample(seq(from = 0, to = 1, by = 0.05), size = 100, replace = TRUE)
# create areas of interest in coarser raster
r2[r2 < 0.9] <- NA

Now i can't seem to find a function (might be blind) that will tell me what i want to know straight off - i thought zonal{raster} might have but it wont return a count function. So this is my work around:

# disaggregate the coarser raster to the same res as the finer raster
r2 <- disaggregate(r2,fact=c(100,100))
# mask the fine by the coarse
r3 <- mask(r1,r2)
# and return a frequency table of the finer resolution
freq(r3,useNA="no")

But this seems a little round the houses.

ISSUE 1: Is there a function to zonal count in R with differing resolution rasters?

ISSUE 2: I up-scaled the above method to 2 very large rasters but got the error "Error: Failure during raster IO", why so?  ## FIXED: issue with PC, not R ##

ISSUE 3: what if I change each cell in the coarser resolution to have an ID instead of a value and I want to count the frequency of values in the finer res raster per ID?

(have also tried changing coarse res raster to binary 1/NA and multiplying but also get the same error as Issue 2 - im using a powerful computer that has worked with bigger data in big stacks, so the issue is not that).


texts: Using spatial join in ArcGIS for Desktop?

texts:
I have two spatial layers of about the same (not exactly same) spatial extent. Layer#1 has about twice the number of polygons than layer#2 (in thousands).

I am trying to calculate values for Field#X and Field#Y present in layer#2 for layer#1 based on spatial location (using right click -> joins and relate, etc.).

What would be the best way to do it?

I am trying to do it using "join data from one layer based on spatial location" and then averaging. But once I sum up the total for Field#X and Field#Y for the resulting layer, it is a much different than the total for the same in the original layer (i.e. layer#2). It should be comparable as the spatial extents are about the same.

Any suggestions on how to address this?


texts:
I have a .tiff and a .tfw file that I want to change into a GeoTIFF. When I try to open the file in ArcMap and ENVI 5.0, it is successful. However, in Google Earth it is not.

Is there anyway to create the GeoTiff in ArcMap or ENVI?


texts: QGIS best way for capturing features without geometry and child features

texts:
I have two (PostGIS) layers in my QGIS (3.18) project and I defined a relation between them.
The first layer PARENT has no geometry.
The second layer CHILD is the child layer and has a geometry.
In the layer form properties of PARENT I made the relation visible, because I want to see all children in the attribute form of the PARENT as a list.
Now I want to capture a new PARENT with some CHILD features:

switch PARENT to edit mode
add record (Ctrl +)
fill the form with PARENTs attributes and press OK.
unfortunately the form is closed now. When I want to capture some CHILD features for my new PARENT I have have to search it in the attribute table - what is not so funny, when the PARENT layer has many features.

Is there a better way to it?
Maybe I can prevent the form from closing automatically?


texts: ArcGIS Feature Class verse Feature Layer

texts: In general is it better (more efficient) to work with Feature Classes (FC) or Feature Layers (a feature layer of a FC). I know when I construct model I always (or generally) work with feature layers which got me to think that maybe layers were more efficient (especially with very very large feature classes.


texts: Create multiple new shapefiles using information contained in a table

texts:
I'm hoping to create a script that creates a new shapefile for each record in a table.  My python experience is limited, and it has been a while since I used it, so I'm having a hard time getting started.  My table has fields for ID, longitude, latitude, distance and bearing.

For each record, I need to:

create new shapefile with a filename that matches the ID field
create a series of points within the newly created shapefile at 1 meter increments along a line that is defined by the longitude/latitude (this is the start point), distance and bearing fields.

Here is my initial code:

import arcpy
table = "C:/folder/project.gdb/transect_table"
fields = ["TRANSECT_ID", "LATITUDE", "LONGITUDE", "Heading", "Distance"]
shp_path = "C:/folder/output_shapefiles"
with arcpy.da.SearchCursor(table, fields) as cursor:
    for row in cursor:
        arcpy.CreateFeatureclass_management(shp_path, cursor[0])

This bit of code is creating the collection of empty shapefiles with filenames being pulled from the ID field of the table.

Can anybody suggest an approach to populating the shapefiles with points?

Based on the input from radouxju, I've updated the code to the following (note that I'm creating feature classes in a gdb as opposed to shapefile; my shapefiles were getting locked when creating the points but feature classes seem to avoid this issue):

import arcpy
import math

path = "C:/folder/Output_FeatureClass.gdb"
table = "C:/folder/Transects.gdb/table"
fields = ["TRANSECT_ID", "LATITUDE", "LONGITUDE", "Heading", "Distance"]
with arcpy.da.SearchCursor(table, fields) as cursor:
    for row in cursor:
        arcpy.CreateFeatureclass_management(path, cursor[0], "Point","","","", "C:/folder/WGS1984.prj")
        with arcpy.da.InsertCursor(path + "/"+  cursor[0], ["SHAPE@XY"]) as insertcursor:
            for i in range(row[4]):
                xy = [(row[2]*i*math.sin(row[3]),row[1]*i*math.cos(row[3])),]
                insertcursor.insertRow(xy)

This is working pretty well, as all of my feature classes are getting created and are populated with the appropriate number of points.  However, the distance between the points is not the desired 1 meter (it's closer to 15,000 km!).  I suspect that I might have an issue with mixing units, as the input positions are in degrees, and the offset between points should be in meters, but I'm not quite sure.

Any suggestions?


texts: How to add more events like Control.Click and Handler.Point for a single variable in OpenLayers

texts:
I have two points which are marked dynamically on the map and when I mark the start point
it is returning lon,lat, but I want to assign this coordinates to the variable that marked/handled the start symbol and a point.
Ex:

measureControls={
              start: new OpenLayers.Control.DrawFeature(start, OpenLayers.Handler.Point),
              stop: new OpenLayers.Control.DrawFeature(stop, OpenLayers.Handler.Point),

for the same start and stop, I would like to add event 'click' and assign the point when clicked.

Thank you


texts:
When I open QGIS, the Layers Panel defaults to a width that doesn't show the "Remove Layer/Group" button. So I have to resize the panel to show the button every time I open QGIS.

Is it possible to set a default width for this panel so I don't have to resize it?

I am using QGIS 2.18.5 on Windows 10.


texts: Creating a new Normal template in ArcCatalog?

texts:
I need to restore my ArcCatalog template to its original default. I am missing the main tab above the toolbar i.e. edit/file/customize etc...I am running ArcCatalog 10.3 on a Windows 10 operating system.

I understand that I need to delete or rename the existing Normal.gxt file in order to generate a new one.  For the life of me I cannot find the template anywhere.

Where is it?


texts: Use Python to find the pathname of "Current" Arcmap

texts:
I am trying to find the path name of the current ArcMap document that the python script is currently running in.

For Example:

"F:\2012\1111\MyArcMap.mxd"

I tired using

sys.argv[0]

as shown here

Finding script location... dynamically

but this gives me the path of the python script not the .mxd

Any thoughts?


texts: Creating cartogram maps in ArcGIS for Desktop?

texts: I have followed the exact steps from http://www.gislounge.com/how-to-make-area-cartogram-maps-in-arcgis/ and downloaded from arcscripts.esri.com. After adding the cartogram tool into the toolbox and I try to invoke it, it says the tool (newly added cartogram tool) is not valid. And when I reopen the same file, the "+" icon in front of the cartogram tool will be missing as if I can't go to the specification pop-up screen. Could this be a compatibility issue (the Arcgis is 10.2.2)? If so, where can I find the right tool to add?


texts: Unable to compute intersection between two polygons using PyQGIS?

texts:
I have two polygons in different layers. I am trying to find the intersect between the two polygons. Below a figure showing my test case. I initially wanted to find the intersect between the green polygon and a bunch of other polygons and wrote a code for that. For all but one polygon this functioned as expected, shown below is the one polygon for which my function fails.

My code:

from qgis.core import *
from PyQt4.QtCore import *

layer1=None
layer2=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() =='Poly':
        layer1 = lyr
    if lyr.name() =='Buffer':
        layer2 = lyr

features1 = layer1.getFeatures()
features2 = layer2.getFeatures()

for feature1 in features1:
    for feature2 in features2:
        geominter=feature1.geometry().intersection(feature2.geometry())

print(geominter.area())
print(geominter.asPolygon())

layer = QgsVectorLayer('Polygon', 'Intersect' , 'memory')
prov = layer.dataProvider()
layer.startEditing()
feat = QgsFeature()
feat.setGeometry(geominter)
prov.addFeatures([feat])
layer.commitChanges()
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])
layer.triggerRepaint()

This code worked for all other polygons I needed to find the intersect for except this one. The output layer contains one element (as expected) but nothing is plotted. The print outputs of this code are:

30.2517066088
[]

It accurately computes the area of the polygon I expect to obtain but the .asPolygon() output is empty, which I expect is why nothing is plotted. What could be the reason for this, and how can I obtain the intersect between these two polygons?

I am using QGIS version 2.18.11 on Windows 8.