texts: Why is an extra field necessary when creating point shapefile from csv files in Python?

texts:
I have a Python script that creates a point shapefile from a csv file within a Python script. This script gives me a ValueError for the longitude field (LON) if there isn't an extra nonsensical field to the right of it. Below is the error and how the csv file has been formatted to remove this error:

    lonValueIndex = valueList.index("LON")
    ValueError: list.index(x): x not in list

DATE       LAT      LON          TEST
2/19/2011  34.27531  -118.21071    a
2/19/2011  34.18069  -118.34079    b
2/19/2011  34.0215   -118.21857    c

Do you see that Test field with values of a, b, and c? If I don't put a field to the right of the LON field I get the ValueError. Below is part of my script. Does anyone know why that extra field is necessary? Thank you.

import arcpy, csv
arcpy.env.overwriteOutput = True

#Set variables
arcpy.env.workspace = "C:\\GIS\\StackEx\\"
outFolder = arcpy.env.workspace
pointFC = "art2.shp"
coordSys = "C:\\Program Files\\ArcGIS\\Desktop10.0\\Coordinate Systems" + \
           "\\Geographic Coordinate Systems\\World\\WGS 1984.prj"
csvFile = "C:\\GIS\\StackEx\\chicken.csv"
fieldName = "DATE1"

#Create shapefile and add field
arcpy.CreateFeatureclass_management(outFolder, pointFC, "POINT", "", "", "", coordSys)
arcpy.AddField_management(pointFC, fieldName, "TEXT","","", 10)

gpsTrack = open(csvFile, "r")

headerLine = gpsTrack.readline()
#print headerLine
valueList = headerLine.split(",")
print valueList
latValueIndex = valueList.index("LAT")
lonValueIndex = valueList.index("LON")
dateValueIndex = valueList.index("DATE")

# Read each line in csv file
cursor = arcpy.InsertCursor(pointFC)
for point in gpsTrack.readlines():

   segmentedPoint = point.split(",")
   # Get the lat/lon values of the current reading
   latValue = segmentedPoint[latValueIndex]
   lonValue = segmentedPoint[lonValueIndex]
   dateValue = segmentedPoint[dateValueIndex]
   vertex = arcpy.CreateObject("Point")
   vertex.X = lonValue
   vertex.Y = latValue
   feature = cursor.newRow()
   feature.shape = vertex
   feature.DATE1 = dateValue
   cursor.insertRow(feature)

del cursor


texts: Using Update Cursor to Calculate Field from List of other Values

texts:
I'm working on a script that will use a predefined list of strings to calculate the value of a new field.  The idea is this: if the value of field 'EF_Type' is equal to the value of a string in the field 'Schools', then the new field 'EF_Group' will be calculated to equal 'Schools.'

    School = ['CollegesUniversities', 'PrivateSchools', 'PublicSchools', 'SupplementalColleges']
    EmergencyServices = ['EMS_Stations', 'Hospitals', 'UrgentCare', 'RedCross']

for Sublist in [PointsInt, LinesInt]:
for fc in Sublist:
    fcname = os.path.basename(str(fc.getOutput(0)))
    with arcpy.da.UpdateCursor(fc, ['EF_Type', 'EF_Group']) as cursor:
        for row in cursor:
            val = str(row[0])
            if val in School:
                row[1] = 'School'
                cursor.updateRow(row)
            elif val in EmergencyServices:
                row[1] ='Emergency Services'
                cursor.updateRow(row)

When I run my script, no error is thrown, but the second field does not update with the new value.  What am I doing wrong?  I've based this off answers to similar questions here, but can't seem to troubleshoot the issue.


texts: Seeking Geospatial CMS like Cartaro without Internet connection?

texts:
I want to install geospatial CMS like Cartaro and I'm interested in whether this CMS works in an intranet environment without taking data from services such as Openlayers, Google map, etc.

Is there a database (world map) as the demo site http://demo.cartaro.org/  that can be used offline?


texts: ArcMap - Adding raster data automatically into geodatabase feature class

texts: I have a point feature class in a file geodatabase, I created a raster column to add some photos, but since I have 1300+ lines I can't add every single photo manually, is there a way to add them automatically ? I have a column with the the photos names.


texts: Calculate Field in python, using codeblock. Indentation Errors

texts:
I have been trying to calculate a field using code block in a python script. For this I am trying to set a variable for the code block called "codeblock".  What I ended up finding was some interesting indentation errors before finally getting it to work. Below are two examples of code, one that didn't work and one that did work. I would like to find out why one works and the other does not.

Example 1:(DID WORK)

Codeblock = """def RemoveNull(X):
    if X is None:
        return 0
    else:
        return X"""

Example 2:(DID NOT WORK)

Codeblock = """
def RemoveNull(X):
    if X is None:
        return 0
    else:
        return X"""

The only difference between the two was where I placed the "def RemoveNull(X):"

Does anyone know why example 1 worked and example 2 did not? I would guess it has to do with the triple quotation marks, but I am unsure.


texts: How to preserve old range of values clipping multi-band raster in QGIS?

texts:
Old raster:

New one after use of "clipper" or GDAL clipping:

Is it any other (faster) method than "save/load style" to preserve previous values of raster? I'm running QGIS 2.14.


texts: Free Web-service providing Satellite Images

texts: Is there any web-service currently available which is hosting satellite images as open-source, just like OSM web-service hosting vector maps? We want to make a real-time traffic monitoring commercial application for Turkey and want a hybrid map as our base-layer (OSM layer + Satellite Images with 50% transparency). Images where roads are easily visible are suffice and young images are more preferable.


texts: Select by Attribute all features with values that are multiples of 5?

texts:
In my shapefile, I have a bunch of features and one field is a integer numerical attribute. I'd like to select all of the features whose value is a multiple of another value. For instance, I have a field with typical values of 100, 101, 102, 103, 104, 105, 106, etc. (up to 1000) and I would like to select all features with a value that is a multiple of 5.

I am using ArcGIS Desktop 10.


texts: Lightweight python modules to convert LAS data formats to text formats?

texts:
I searched on the net and found pyLAS, but it seems to require quite a few other modules and I couldn't get the install right.

Any pointers will be much appreciated.


texts: Error in proj4js-Projection while upgrading from OpenLayers 5 to OpenLayers 6 (ETRS:25832 projection)?

texts:
I have a map with some layers in Projection EPSG:25832. This works well with OpenLayers 5.3.3. But when I change to OpenLayers 6.4.3 (or all others down to 6.0.1) then there is an error in the Definition of my projection. I am using proj4js in the version 2.6.2 (also tried older versions).

https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.6.2/proj4.js

Error thrown from Proj4js:
Uncaught TypeError: coordinates must be finite numbers .... proj4-2.6.2.js:1:4846
Sourcecode in Openlayers-5 and Openlayers-6
var oBoundsGesamt = [ 456700,5200500,773480,5835251 ];

proj4.defs(sEpsgNummer, "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs ");
ol.proj.proj4.register(proj4);

var oOlProjection = ol.proj.get('EPSG:25832');
oOlProjection.setExtent( oBoundsGesamt  );

var oMapObj = new ol.Map({
   target: 'ebMapObj',
   layers: [ ... ],
   view: new ol.View({
        projection : oOlProjection,
        extent: oBoundsGesamt,

Can somebody help me with this? The error is thrown at the end when I try to zoom to an Extent:

ol.control.ZoomToExtent({ extent: oBoundsGesamt }) or
oMapObj.getView().fit(oBoundsStart, oMapObj.getSize());


texts: Use Common Table Expressions with pgr_TSP

texts:
I have a query that uses a handful of CTEs to return a generated set of points. The actual query is more complicated, but the following will do as an example:
WITH geom AS (
    SELECT unnest('{
        POINT(-113.3997627 53.5721641),
        POINT(-113.4037653 53.452761),
        POINT(-113.7066278 53.4571305),
        POINT(-113.7860125 53.5788978)
    }'::TEXT[]) AS point
    ORDER BY 1
),
points AS (
    SELECT (ROW_NUMBER() OVER ())::int AS id, point::geometry AS point
    FROM geom
)
SELECT *
FROM points;

I would like to use these points in the travelling salesperson problem function, but I haven't been able to make a simple query that achieves this.
Using the pgr_tsp function at the end of the query fails with at relation "points" does not exist error:
WITH geom AS (
    SELECT unnest('{
        POINT(-113.3997627 53.5721641),
        POINT(-113.4037653 53.452761),
        POINT(-113.7066278 53.4571305),
        POINT(-113.7860125 53.5788978)
    }'::TEXT[]) AS point
    ORDER BY 1
),
points AS (
    SELECT (ROW_NUMBER() OVER ())::int AS id, point::geometry AS point
    FROM geom
)
SELECT seq, id1, id2, round(cost::numeric, 2) AS cost
FROM pgr_tsp('SELECT id, ST_X(point) AS x, ST_Y(point) AS y FROM points', 1);

Moving the entire query into the first parameter of the function seems to work, but I would need to either duplicate the CTE expressions outside of the function or store them in a temporary table to join them in the resulting order:
SELECT seq, id1, id2, round(cost::numeric, 2) AS cost
FROM pgr_tsp($$
    WITH geom AS (
        SELECT unnest('{
            POINT(-113.3997627 53.5721641),
            POINT(-113.4037653 53.452761),
            POINT(-113.7066278 53.4571305),
            POINT(-113.7860125 53.5788978)
        }'::TEXT[]) AS point
        ORDER BY 1
    ),
    points AS (
        SELECT (ROW_NUMBER() OVER ())::int AS id, point::geometry AS point
        FROM geom
    )
    SELECT id, ST_X(point) AS x, ST_Y(point) AS y
    FROM points
$$, 1);

Is there any workaround for what I am trying to achieve, or should I use a transaction with temporary table?
I am using psql9.5, postgis2.2.1, and pgrouting 2.1.0. Yes, I am planning on updating them.


texts: Problem with ST_IsValid: Unknown geometry type: 9 - CompoundCurve

texts:
The following is an example of a gml:Curve that was loaded into Postgres using FME.

<gml:Curve srsName="urn:ogc:def:crs:EPSG::27700" gml:id="os.vmd.ny.38124.geom">
<gml:segments>
<gml:LineStringSegment>
<gml:posList>379352.97 580000 379387.02 580022.77 379392.7 580028.98 379412 580031 379419.11 580023.97 379420 580016.95 379446.92 580000</gml:posList>
</gml:LineStringSegment>
<gml:LineStringSegment>
<gml:posList>379446.92 580000 379508 580004.5 379521 580016.5 379529.14 580047.36 379531.91 580057.92 379552.03 580087.79 379559.12 580095.84 379606.5 580115.5 379618 580112.5 379630.79 580097.79 379636.51 580091.2 379667.5 580050 379717.24 580096.52 379729.15 580088.65</gml:posList>
</gml:LineStringSegment>
</gml:segments>
</gml:Curve>

PostgresSQL accepted the geometry but all attempts to do anything with it yielded the error:

  ERROR:  Unknown geometry type: 9 - CompoundCurve

including ST_IsValid(geom). Therefore, there is no way to even find out the unique id of the offending geometry, as

SELECT id FROM water WHERE ST_IsValid(geom) = 'f'

returns this error rather than the id.

Is there anyway around this?

Update following comments: It should be noted that all functions return this error, including ST_GeometryType(geom), and ST_AsText(geom). @Vince suggested using a binary search to find the offending rows, which would work, but with millions of rows and hundreds of invalid geometries, would not be fast.

As a related question is there any reason why ST_IsValid returns an error, rather than regarding an error as valid = false.

Fortunately, I have a shapefile of the same area, in which the geometry is encoded as a MultiLinestring, but I would like to avoid having to reload the entire table if possible.


texts: Looking for an android solution for offline kml/kmz file viewer with gps support

texts:
I have some KMZ files of cadastral maps. I need to fiscally go to the boundaries of my lands (parcels) and take a look at where they are with the gps and maps support. Some of these places are not covered by the phone network signal.

For this reason I'm looking for an android app that will import my kmz files and let me reach the boundaries position with gps help.

I've used google earth app, but gps position is not working.

Google maps engine app, instead, works well but only with online maps.

Any idea about an app to work offline?


texts: Importing Data Dictionary (DDF) into QGIS (Fields)

texts:
Are there any tips/tricks on getting all the Attributes from a data dictionary for a GPS unit to sync with the Fields in a QGIS style?

Right now, in order to update the drop-down options in QGIS, I have had to semi-manually enter and update the Fields and Values when any of the data dictionary items change (Trimble). The default "Widget" type for imported shapefiles in QGIS is "Text Edit," I have to change that to "Value Map" or other widget type and import the fields/values via csv for each field. Any way to automate that process by importing a DDF or other method?  I have ~100 Attributes/Field types.


texts: Collecting simultaneous features or reducing tracking interval on location tracking in ArcGIS Collector?

texts:
I encountered a problem in using collector for our data collection on the field. We are archaeologists doing prospection work and needing to create point features for found artefacts as well a streaming line features (with somewhere between 1-5 seconds tracking intervals) to know exactly where we walk in surveyed zones. We need this precision in our line features to calculate prospection intensity in said surveyed zones.

But, Collector seems to be unable to stream a line feature while you are taking points. The solution a fell upon was the use of the esri layer named « location tracking » which leaves a bread crumbs trail of points where you go with an interval of 30 seconds at the least.

Is there a way to modify the ESRI location tracking layer for it to use an interval under 5 seconds?

Or, by any chance, did anybody found a way to simply stream a classic line feature while collecting points?


texts: Extracting time series from image collection in Google Earth Engine?

texts:
I tried to perform a pixel based time series change detection algorithm in GEE for a defined region, so i intended to extract every pixel time series values into something like a table or matrix, and implement the pixel based algorithm by rows or columns.

I find out a similar question and rewrite the code by following the solution given by Nicholas Clinton as follows:

#define the interested path and row to select the images
path = 120
row = 32

#define the time period to filter the images you want
start = ee.Date.fromYMD(1983, 1, 1)
finish = ee.Date.fromYMD(2019, 1, 1)

#rename the OLI images
l8_bandlist = ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'pixel_qa'])
rename_list = ee.List(['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'tbb', 'qa'])

#build the image collection
l8_sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filter(\
ee.Filter.eq('WRS_PATH', path)).filter(ee.Filter.eq('WRS_ROW', row)).filterDate(\
start, finish).select(l8_bandlist, rename_list).sort('system:time_start')

# using shapefile to restrict the area of interest
shp = ee.FeatureCollection('users/myzhenghrsc/testarea').geometry()

#reduce the region restricted by uploaded shapfile into points correspond to OLI pixels
dictionary = ee.Image.pixelLonLat().reduceRegion(
reducer=ee.Reducer.toCollection(['longitude', 'latitude']),
geometry=shp,
scale=30)

#function to build the piont featureCollection
def __rebuildpoints(feature):
lon = feature.get('longitude'),
lat = feature.get('latitude'),
return ee.Feature(ee.Geometry.Point([lon, lat])), {
    'featureID':
    ee.Number(lon).multiply(1000).round().format('%5.0f').cat('_').cat(
        ee.Number(lat).multiply(1000).round().format('%5.0f'))
}

#perform the function above
points = ee.FeatureCollection(dictionary.get('features')).map(__rebuildpoints)

def triplets(image):
  def feature(value):
    return value.set({
        'imageID': pointvalue.id(),
        'timeMillis': pointvalue.get('system:time_start')
    })

  pointvalue = image.reduceRegion({
    'collection': points,
    'reducer': ee.Reducer.first().setOutputs(image.bandNames()),
    'scale': 30
})
return pointvalue

turplets = l8_sr.select('red').map(triplets)

def newformat(table, rowId, colId, rowProperty, colProperty):
rows = table.distinct(rowId),
joined = ee.Join.saveAll('matches').apply({
    primary=rows,
    secondary=table,
    condition=ee.Filter.equals({
        leftField=rowId,
        rightField=colId
    })
})

def __row(row):
    def __feature(feature):
        feature = ee.Feature(feature),
        return [feature.get(colId), feature.get(colProperty)].flatten()

    values = ee.List(row.get('matches')).map(__feature)
    return row.select([rowId, rowProperty]).set(ee.Dictionary(values))

return joined.map(__row)

results = newformat(triplets, 'imageID', 'featureID', 'timemillis', 'RED')

But i got a error as follows when i perform this code:

  Invalid argument specified for ee.Number(): ee.ComputedObject({
    "type": "Invocation",
    "arguments": {
      "object": {
        "type": "ArgumentRef",
        "value": null
      },
      "property": "longitude"
    },
    "functionName": "Element.get"
  })


texts: In ArcGIS 9.3, how do I eliminate the slivers of missing data between polygons?

texts:
I am editing a single-band .tiff file in ArcGIS 9.3.  I started with an image that I had to georeference, but unfortunately, the image had these little slivers of missing data between the data cells.  I get left with a raster file that has lots of "slivers" of nodata between my coding cells (shown here in white).

I would like to use something analogous to vector snapping, in which I eliminate those white slivers, instead making those sliver raster cells equal to the value of one of the adjacent cells.  thank you.


texts: Geometry in PostGIS Table wrong SRID and type

texts:
I have loaded some points into a table as Multipoint, they have come through with an SRID of 0. I need to change these to have a SRID of 7844 and also change them to be of type Point. There are other records already in the table with type of Point and SRID of 7844.

I have seen how to change the SRID from unknown to another projection here, https://postgis.net/2013/08/30/tip_ST_Set_or_Transform/ using ST_Set.

I have also seen that in PostGIS you can use ST_Dump for multipart geometries. I have also just found Multipoint vs point PostGIS

How do I change specific records SRID and downgrade from multi to individual points in PostGIS in SQL without affecting the other records.


texts: Error in points or military academy dataset?

texts:
I travel regularly to sites around Vietnam reported in declassified Australian military documents. All grid coordinates are reported in abbreviated MGRS and are usually quite accurate.

In recent years a university has created a large resource intensive map of most of the same sites which appear on Google Maps different positions. The explanation I had from the uni was that the war era coordinates were using Indian 1960 datum and they've had to correct this to fit into WGS.

However, back in those days they were not using GPS, but reading from topo maps with a MGRS grid overlay. I'm no expert GIS wiz, just an end user, but I would think that the MGRS grid would be the same as WGS / UTM and the grid coordinates presented in the declassified documents would not need transformation.

In all my years of using the grids to overlay georeferenced sketch maps onto Google Earth etc the landscape features have always matched up.

Below I present one example of a single point reference.

Can anyone explain the discrepancy with the university map?

I'd hate to think I'd been out of calibration all these years.

The original quoted location is YS508549. Naturally this is an abbreviated MGRS.

If I expand this I believe it should look like 48PYS5080054900;

Or 750800.00 m E 1154900.00 m N (10.439367°N 107.290986°E decimal degrees)

When I plot this onto Google Earth it appears just south of Duong Bach Mai school, some 4-500m away from the location on the university battlemap (diagonally opposite side of the main hwy 44 /Vo Thi Sau road)

https://vietnam.unsw.adfa.edu.au/battlemap/?contact-filter=contact-units-involved:6866,6950,7011,7019,8629,8723&basemap=mapbox-terrain&layers=contact-individual&timeline=1969-05-26,1970-04-17&incident=2737&at=10.456847037511139,107.29440676845547,14


texts: Getting points with lowest value and exact location for each cell?

texts:
I have an assignment to create a simplified depth scheme.
The calculation has to be done for a large area automatically.

There are depth points every 0.5 meters (marked black in the picture). By using “point to raster” function I can generate a raster layer with defined size cells (e.g. 4x4 m) and also request to calculate lowest (shallowest) values for each cell.

Then I can create points layer from the raster generated before, which results is a layer of depth every 4 m (points in red) with lowest values.
However the points are generated in the center of the cell. They do not represent the actual place within the cell with lowest value (yellow points).

How can I get the points with lowest value and exact location for each cell?


texts:
Possible Duplicate:
  Comprehensive list of Latitude and Longitude values of cities round the world

i would like to know if is out there any downloadable database of world locations,

i know there is openstreetmap db available, but i would like somenthing lighter then that, i just would like a world cities + region + nation + coordinates db, cause my intent is just to query the db to find near places to a given city.

So doesn't matter if db has cities with 1K population for example, i don't need them all, i just need the more/main important cities up to 20/30K citizens i mean and their related coordinates + region + nation , thanks i know it should be like it there isn't out there :P

Maybe someone could help me then on querying openstreetmap db to remove all not needed data !?


texts: How to export multiple class training data from ArcGis into R?

texts:
I want to perform some machine learning classification of images using R. I'm having issues finding out how to export my training data from ArcGis into R and I'm looking for some help. Currently I do not have access to ArcGis but I will provide visual examples of what I am trying to achieve.

I have an image (with multiple landsat layers) and some points already classified, as the following image shows.

I want to select polygons around that image that I know are the same class as the points, like in the following image:

And then I want to export these classes as training data for R, so I can classify the rest of the image. How can I achieve this? I've been searching in the gis stackexchange site and I found a similar answer (in this case for random forest classification) in https://gis.stackexchange.com/a/39103/48810, my issue is regarding the part

# READ POINT SHAPEFILE TRAINING DATA
sdata <- readOGR(dsn=getwd() layer="inshape")

where only one shapefile is read? I have also found another answer in https://gis.stackexchange.com/a/26916/48810 where multiple shapefiles are apparently exported, but I don't know if it's possible to import multiple shapefiles into R. So I'm asking for help on the proper way to do this kind of export and classification using ArcGis and R.

The arcGis version I will be using is the 10.2.2 one, in case it matters. If necessary I will provide additional information.


texts: Access secure services from ArcGIS Portal with API JavaScript

texts: I was reading this tutorial, that explains how to access secured services on ArcGIS Online. But I wonder if there is a way to do the same but with Portal and access my secured own data.


texts:
I have geopandas with CRS in meters (ITM {'init' :'epsg:2039'}).

I converted the geometry from espg:4326 so I have very high precision which I don't require: POINT (177559.4776412862 594546.4983892406)

How can I truncate the redundant digits after the floating-point?
I tried to use:

df.geometry = df.geometry.simplify(0.5) but I get an error:

  can't set attribute


texts: Removing ONLY black background in QGIS (without removing black dots within the image)

texts:
Related to question: Removing black background and black font from raster in QGIS

I have also a black backround and when doing as it is suggested in the answer the black backround goes away but also points from the image within the region I want to keep are also erased. There may be black dots which dissapear when doing this. So white dots appear.  I thought the option "clip the extent of the target dataset to the extent of the cutline" should control from which line on the nodata should be applied but it isn´t the case. I have tried all the options (0, 100, leaving it empty, ...). Is there another thing I could do?


texts:
I have a virtual environment, (my-virtualenv) bash-3.2$, running with python 3.6.5 set up using pyenv.

I installed the modules listed on the readme (owslib, PyYaml, psycopg2, jinja2, pygments, numpy, plotly). My system has python 2.7 and python 3.7 outside the virtual environment.

How do I install QGIS just within the virt environment?


texts: Missing geometry in Sentinel 2 footprints

texts:
I've noticed some scenes on bucket that contain empty (no geometry) detector footprint masks. Examples:
http://sentinel-s2-l1c.s3.amazonaws.com/tiles/11/S/MT/2016/8/17/0/qi/MSK_DETFOO_B03.gml
http://sentinel-s2-l1c.s3.amazonaws.com/tiles/36/L/XR/2017/11/10/0/qi/MSK_DETFOO_B04.gml

I didn't find it mentioned in Sentinel 2 quality report. I'm sorry to bother, I know it's not the right place to report an issue. Where should I report issues related to Sentinel 2 on AWS?

UPDATE:
I have checked the masks at Copernicus SciHub and they are empty as well. So I wrote to EO Support.


texts: What is the height reference point when using Longitude, Latitude , Height coordinates?

texts:
I'm using a library to convert LLH(WGS84) to ECEF coordinates.   Does height (for corresponding longitude and latitude) denote height above/below the designated earth radius (I'm using 6378137.0m) or distance from center of the earth?
Edit Code Refs:
I'm basing my code off this repo
But converting to C++.  However I don't believe that code is very accurate.
Looking at code snippets behind this web page
Exploring pymap3d but it's only for Python 3.x and I require 2.x
Edit more info:
Also from reading Geographic coordinate system it states:

A common choice of coordinates is latitude, longitude and elevation

And then defines Elevation as:

The elevation of a geographic location is its height above or below a
fixed reference point, most commonly a reference geoid, a mathematical
model of the Earth's sea level

I'm attempting to convert from ECEF to ENU coordinates using the above code/libraries.  ECEF to ENU requires a reference LLH but can't find any documentation as to what height value is expected, distance from center of earth or elevation above sea level.
So to clarify my question, if one communicates a location using Latitude and Longitude, is the "height" standardized as elevation above sea level?


texts: How to make projected geometry and update it to postgresql?

texts:
I have an excel file containing some unprojected geographical coordinate values such as:

Time,Lat,Lon,Alt
2017:11:21:13:44:07:132,61.8984838,24.2221901,165

I wrote a java code to store the data to postgresql database. Now, my postgres database contains a table named 'info' which contains all the data from excel file. I also added a gid column to give unique serial number for all data and a geometry column to my database with null value. I want to put projected UTM value to that geometry column. My required UTM zone is WGS 84 UTM 35N and its EPSG code is 32635. Is st_transform able to do this projection? My sample query is here:

Update info set geom = sq.geom
from (
    select ST_Transform(geom,32635) as geom
    from (
        SELECT ST_setsrid(geom_v,4326) as geom
        from (
            select ST_MakePoint(lon,lat,alt) as geom_v
            from (
                select latitude as lat, longitude as lon, altitude as alt from info order by gid asc
            ) as fpp
        )as ftt
    ) as tr
) as sq

I added st_setsrid as without this, postgres was returning "undefined 0 srid values" error. I set the srid to 4326 just for an example. It could be any valid srid. Anyway the query is filling the geom column with same value for each row. Such as:

geom
geometry
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F000026218A29179C154199C0063B12315A41'

But if I try this:

select ST_Transform(geom,32635) as geom
from (
        SELECT ST_setsrid(geom_v,4326) as geom
        from (
            select ST_MakePoint(lon,lat,alt) as geom_v
            from (
                select latitude as lat, longitude as lon, altitude as alt from info order by gid asc limit 10
            ) as fpp
        )as ftt
) as tr

Then, it gives 10 different geometry values. Such as:

geom
geometry
'01010000207B7F000026218A29179C154199C0063B12315A41'
'01010000207B7F0000463A245D169C1541B86D933712315A41'
'01010000207B7F0000629183D5149C1541EBB6462A12315A41'
'01010000207B7F00007ECE9CA0139C15418EC6E61712315A41'
'01010000207B7F000064C89D65129C15415AD97AF911315A41'
'01010000207B7F0000B0693E76119C1541291311E011315A41'
'01010000207B7F0000F868ED7D109C154120819BC411315A41'
'01010000207B7F00003A478EDD0F9C1541E15E44AD11315A41'
'01010000207B7F00003A478EDD0F9C1541E15E44AD11315A41'
'01010000207B7F00003A478EDD0F9C1541E15E44AD11315A41'

Why this?


texts:
What are some ways to embed a Leaflet map onto one's website?

I've explored other posts that mention iframe, and plugins like MapsMarker.

But, what is a step-by-step way to do it?

My site uses a Wordpress CMS. Can I host my html/css/js files directly on WP? Where do I put them? Or, do I need to make a folder on the file manager of my host (Host Gator)?

(I've been self-hosting in MAMP to test it out, which works well).


texts: Is it possible / wise to work only with .TAB files in QGIS

texts:
QGIS is great for it's ability to open and project almost anything I throw at it. I'm thinking of moving from Mapinfo 10.5 to QGIS, is it possible to just keep using the same .TAB files or would this limit certain features of QGIS or cause other problems?

The main reason I'd like to keep using the .TAB files is that other members of the office are sticking with MapInfo and I don't want to cause an office data schism.


texts: How do I classify number values in Python to assign it to a group?

texts:
I have calculated the number of points that exists for a SqMi area and am trying now to assign it to a certain group based on where it falls in the classification. Does anyone know a tool that I can run that will generate an output, or the python code that will help me to test these values to see where they need to be. I have inserted a sample code and I will write some pseudo code so that you can follow my logic to help.

#Calculate per sq mi the number of incidents within the zone (Top=15+, High = 12-15, Some = #6-12, Low = >6)
# Get the shape_area number from the CurrentZoneLayer
field = "SHAPE_Area"
with arcpy.da.SearchCursor("CurrentZoneLayer",field)as cursor:
    for row in cursor:
        zoneSqMi = (row[0]/2589988.11)
        print "Zone SqMi: "+ str(zoneSqMi)

# Divide the zoneSqMi by the number selected to get points by SqMi
pointsPerSqMi = numSelected/zoneSqMi
print "Points Per Sq Mi: " + str(pointsPerSqMi)

# Create a series of tests to see which group the points will fall under Top, High,  #Some, Low
# if pointsPerSqMi > 15 then print "Top"
# if pointsPerSqMi <15 and pointsPerSqMi

Eventually these categories will be written into a field value as well, but for now I am just trying to figure out the logic on how to assign it to a group within a script.


texts: Does anyone know where to download a shapefile containing the locations of US interstate and highway exits?


texts: QGIS Grid Division/Grid Product etc. results are wrong projection

texts:
I have several raster grids, all 4320 columns, 2160 rows in WGS 84 projection. They are loaded from ESRI ascii files.
When I use the Processing Toolbox to carry out arithmetical transformations on the grids (e.g. multiply one by another, divide one by another) the resulting grid is correct, but offset from the original projection by some number of pixels. The layer purports to be in WGS 84, but clearly isn't.
Currently I am working round the problem by using the raster package in R, but it seems to me that the QGIS functions are broken (unless I am missing something obvious).
Any suggestions gratefully received.


texts: Getting parameter menu in QGIS modeler to appear?

texts:
I lost, by accident, the parameter menu in QGIS 2.14 Graphic Modeler.

Is there any tip for how can i get it back.


texts: QGIS crashes frequently since upgrading to 2.18.2 on Mac

texts:
I upgraded QGIS to 2.18.2 from here http://www.kyngchaos.com/software/qgis.

I did note this warning there:

WARNING: QGIS will crash if Qt4 developer components are installed in the standard /Developer location. Either rename /Developer/Applications/Qt/Plugins or uninstall Qt before running QGIS.

But I don't see Qt4 installed there.  Maybe I'm not searching correctly, as I'm still somewhat new to the Mac paradigm and finding the install locations of software.

The crashes happen most often when I'm dragging or zooming a map, and I get a Mac problem report that includes this information (I'm not pasting the entire report because it is pretty long):

I'm curious if others have encountered this and how you resolved it?

Process:               QGIS [44719]
Path:                  /Applications/QGIS.app/Contents/MacOS/QGIS
Identifier:            org.qgis.qgis2
Version:               2.18.2 (2.18.2 [exported])
Code Type:             X86-64 (Native)
Parent Process:        ??? [1]
Responsible:           QGIS [44719]
User ID:               501

Date/Time:             2017-03-02 09:16:43.972 -0800
OS Version:            Mac OS X 10.10.5 (14F27)
Report Version:        11
Anonymous UUID:        4C33FF15-D5BB-A25F-93B6-A78DE682F1A7

Sleep/Wake UUID:       6F832B34-2C1A-48B1-88B5-58CC530EFAE8

Time Awake Since Boot: 1700000 seconds
Time Since Wake:       6500 seconds

Crashed Thread:        18  Thread (pooled)

Exception Type:        EXC_BAD_ACCESS (SIGSEGV)
Exception Codes:       EXC_I386_GPFLT

Thread 0:: Dispatch queue: com.apple.main-thread
0   QtGui                           0x000000010fc0eb58 QIcon::~QIcon() + 24
1   QtGui                           0x000000010ff669fe QDockWidgetTitleButton::sizeHint() const + 238
2   QtGui                           0x000000010ff66c12 QDockWidgetLayout::titleHeight() const + 194
3   QtGui                           0x000000010ff67678 QDockWidgetLayout::sizeFromContent(QSize const&, bool) const + 856
4   QtGui                           0x000000010ff6aa82 QDockWidgetItem::minimumSize() const + 146
5   QtGui                           0x000000010ff6e9b1 QDockAreaLayoutItem::minimumSize() const + 97
6   QtGui                           0x000000010ff6e506 QDockAreaLayoutInfo::minimumSize() const + 166
7   QtGui                           0x000000010ff6ea12 QDockAreaLayoutItem::minimumSize() const + 194
8   QtGui                           0x000000010ff6e506 QDockAreaLayoutInfo::minimumSize() const + 166
9   QtGui                           0x000000010ff6e738 QDockAreaLayout::minimumSize() const + 200
10  QtGui                           0x000000010ff9b6e6 QMainWindowLayoutState::minimumSize() const + 38


texts: How to get the nearest points from pgrouting linestring result?

texts:
I'm trying to get the points or landmarks nearest to the result of pgrouting linestring (the shortest path). I'm using OSM map entities and working with planet_osm_point. Here's my code:

SELECT st_asgeojson( ST_Collect( nearRank.nearGeom ) ) As geojson
 FROM ( SELECT ST_Transform( marks.way, 4326 )  as nearGeom, name,
              ST_Distance( road, marks.way )/1000.0 AS distance_km
    FROM planet_osm_point AS marks,
         st_setsrid( ST_GeomFromText('LINESTRING(121.1527761
                     13.9413653,121.154058 13.9413113,121.1550161
                     13.9412372,121.1559079 13.9411891)') , 900913 ) as road
    WHERE marks.amenity <> '' AND marks.name <> ''
    ORDER BY  ST_Line_Locate_Point(road, marks.way),
    marks.way <-> road LIMIT 15 ) As nearRank

The above sql script is not giving the right nearest nodes along the linestring. The linestring above is the result of pgrouting functions. Is there anyone that can tell me what I have missed? The output should be points representing the landmarks near the shortest path.


texts:
I am working on a toolset, and one of the inputs is a field of a shapefile. I have defined the field and I have used it in multiple functions, and it is working without any error.
The problem is that when I want to use the field in an expression in CalculateField, it does not recognize the defined field, but when I use the name of field (name in attribute table) it works.
I think that I make mistake in writing the code for the expression part.

import os, sys
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

###### Define variables
shapes = sys.argv[1]
outputfolder = sys.argv[2]
outputname = sys.argv[3]
field = sys.argv[4]
output = os.path.join(outputfolder,outputname)
arcpy.AddMessage(shapes)
arcpy.AddMessage(field)
arcpy.AddMessage(output)

shapes1=arcpy.AddField_management(shapes,'Classific','TEXT')
shapes2=arcpy.CalculateField_management(shapes1,'Classific',"!Classific!.replace(!Classific!,field)", 'PYTHON3')

Unionfield = arcpy.CalculateField_management(Unionfield, field, "field+ !FIRST_Clas!", 'PYTHON3')
Unionfield = arcpy.CalculateField_management(Unionfield, field, "!Class!+ !FIRST_Clas!", 'PYTHON3')

Also, I have another question about adding the output of my toolset to the contents after each running.
When I run my toolset every time, the output does not appear in the contents, every time I have to import the output from outputfolder.


texts: Make Query Table error: "The value is not a Raster Layer"?

texts:
I'm pretty new in coding and I'm struggling with my MakeQueryTable_managment tool using arcpy. I'm trying to do a query Layer based on a Layer and a table I have both in my gdb.
I keep getting this message error no matter what I change :

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000793: Invalid data element type
ERROR 000840: The value is not a Raster Layer.
Failed to execute (MakeQueryTable)

I don't understand why it's asking for a Raster Layer, as it can be a Table View too as far as I know. When I use the tool with the same Feature class and table in arcMap, it works without any problem and the Layer appears.
Here is the piece of code failing:

MyGDB = arcpy.env.workspace = mygdb.gdb
table_name = "myTable"
myLayer = str(MyGDB) + "myLayer"

arcpy.MakeQueryTable_management([myLayer,table_name],"QueryTable",\
"ADD_VIRTUAL_KEY_FIELD", "", "","myLayer.ID = " + table_name + ".ID" )

Can anyone give me a little trick to solve this?


texts:
I have gone through the following tutorial to start a map server on Fedora 28:

https://wiki.openstreetmap.org/wiki/Mapnik/Installation_on_Fedora_18

Although the tutorial is a bit old, most of the steps work fine. Basically we do the following:

Install and start postgresql-server and postgis
Create and setup database
named gis and user name apache with all privileges granted on geometry_columns and spatial_ref_sys
Build osm2pgsql from source
Get a region data (in osm format)
Import the osm data into gis using ./osm2pgsql -S ./default.style ../region.osm
Install mapnik:
yum install mapnik mapnik-devel mapnik-python proj-epsg
Checkout mapnik data from source:
svn co http://svn.openstreetmap.org/applications/rendering/mapnik
Download coatline data and some other data:

  wget http://tile.openstreetmap.org/world_boundaries-spherical.tgz wget
  http://tile.openstreetmap.org/processed_p.tar.bz2 wget
  http://tile.openstreetmap.org/shoreline_300.tar.bz2 wget
  http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_populated_places.zip
  wget
  http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_boundary_lines_land.zip

Unpack files and setup mapnik environment.

  cd ~/osm/mapnik tar xvzf world_boundaries-spherical.tgz tar xvjf
  processed_p.tar.bz2 -C world_boundaries tar xvjf shoreline_300.tar.bz2
  -C world_boundaries unzip ne_10m-populated-places.zip -d world_boundaries unzip ne_110m-admin-0-boundary-lines_land.zip -d
  world_boundaries cd ~/osm/mapnik cp archive/set-mapnik-env . vim
  set-mapnik-env change "export MAPNIK_DBNAME='osm'" to 'gis' change
  "export MAPNIK_DBUSER=whoami" to 'apache'

Run the script generate_tiles.py (A slightly modified version of it in fact).
Almost all steps are successful on Fedora 28, except the last step which gives me the following error:

  [postgres@localhost mapnik]$ ./z0_generate_tiles.py
  /usr/lib64/python2.7/site-packages/mapnik2/init.py:27:
  DeprecationWarning:  mapnik2 module has been deprecated,
          please use 'import mapnik'    warnings.warn(msg, DeprecationWarning) render_tiles( (-180.0, -90.0, 180.0, 90.0)
  /var/lib/pgsql/mapnik/osm.xml /var/lib/pgsql/mapnik/tiles/ 0 0 World )
  Traceback (most recent call last):   File "./z0_generate_tiles.py",
  line 215, in
      render_tiles(bbox, mapfile, tile_dir, 0, 0, "World")   File "./z0_generate_tiles.py", line 136, in render_tiles
      renderer = RenderThread(tile_dir, mapfile, queue, printLock, maxZoom)   File "./z0_generate_tiles.py", line 64, in init
      mapnik.load_map(self.m, mapfile, True) RuntimeError: XML document not well formed: Entity 'prefix' not defined at line 4064 of
  '/var/lib/pgsql/mapnik/osm.xml'


texts: GeoWebCache WMTS GetTile Request with HTTP POST KVP encoding

texts:
WMTS specifies GetTile requests using HTTP POST, by simply encoding all request params in the body. Is this supported by GWC, and if so, what headers etc. do I have to set?

I naively tried this:

POST http://myserver.com/services/wmts

service=WMTS
request=GetTile
version=1.0.0
layer=mylayer
format=image%2Fpng
tilecol=12
tilerow=34

Which gives me:

<?xml version="1.0" encoding="UTF-8"?>
<ExceptionReport version="1.1.0" xmlns="http://www.opengis.net/ows/1.1"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://geowebcache.org/schema/ows/1.1.0/owsExceptionReport.xsd">
  <Exception exceptionCode="MissingParameterValue" locator="request">
    <ExceptionText>Missing Request parameter</ExceptionText>
  </Exception>
</ExceptionReport>


texts: Limiting memory that raster package uses in R?

texts:
What's the easiest way to limit the memory the raster package uses in R to 5GB of RAM?

Ideally I want any raster that exceeds this size to be written to the hard drive and analysed from a temporary file.

I'm using Windows 7 (64Bit) and R version 3.0.3.

Here's some code I've been playing with on a machine with 16GB of RAM but it keeps maxing out as i'm dealing with large rasters:

require(raster)

memory.limit(size = 5000)
rasterOptions(format="GTiff",overwrite=TRUE,datatype="INT1S",
              tmpdir="C:/Research/BIN",tmptime=1.1,progress="text",chunksize=1000,
              maxmemory=1000)
rasterTmpFile("delete_me_")

r <- raster("myraster.tif")
r[r==0] <- NA

It is at this stage in the script that the RAM usage increases to take all 16GB and then stops working and cancels the function.

What do I need to change here?


texts: Load and convert .txt file to .qml for raster color map

texts:
I have a "QGIS Generated Color Map Export File" in .txt format.  I'm trying to load it using QGIS 2.18 but the only file extension QGIS will accept is .qml file format.  This question relates what I'm trying to do but provides no explanation as to how to load a non-qml file.

The file I'm trying to utilize is located here: http://www.warrenpinnacle.com/prof/SLAMM/GCPLCC/slamm_colors_qgis.txt


texts: How can I make a polygon from intersecting lines of different layers in QGIS?

texts:
I am constructing a geological map which involves constructing boundaries to areas which vary in style across their length depending on how well defined they. To do this I have created two distinct layers called "certain" and "inferred", one solid the other dashed. I wish to automatically create a polygon fill surrounded by the boundaries, but when I use the lines -> polygon tool, I am only given the option to select one layer at a time.

Can I join the two layers together into one but retaining the changing style?

Alternatively, can I possibly change the style of the line within one layer, or can each layer only have one style?


texts: Profile/Transect tool in QGIS like Arc's interpolate transect tool?

texts:
I am quite new to GIS.
I am trying to determine the values along a transect of a raster in QGIS. I can perform this task in ArcMap using the profile tool however I've tried several methods in QGIS (profile tool, point sampling tool) but do not seem to work!
FYI, I'm looking at LST data and using QGIS 2.18.15 for Mac.


texts: RGB TIFF goes black when zoomed out in QGIS?

texts:
Am importing a RGB TIFF into QGIS 3.4 and when I zoom in I can see the image properly, however, when I zoom out the layer turns black (most of it) and red (a small portion). After doing some search I found this, which suggests that this problem is caused by issues with pyramids and statistics with the layer. There is some support there but that is for ArcGIS. I tried building the pyramids in QGIS Raster-->Miscellaneous-->Build Overviews (Pyramids) (I have done it in the past for .asc files), but got the following error

ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.
ERROR 1: An error occurred while writing a dirty block from GDALRasterBand::RasterIO

Does anyone know how I could to build pyramids for a TIFF in QGIS?

If so, what "Overview levels" should I select?

Or does anyone know how to solve this issue for my RGB TIFF?


texts: How to install and configure QGIS WEB CLIENT in Windows

texts:
I have the problem configuring my QGIS web Client. I'm new in QGIS. These are the tool I use:

Xampp
QGIS Server
QGIS Web Client
QGIS Desktop


texts: Python: find a method to calculate the “inner centroid” (also known as labelPoint in arcmap) of a polygon


texts: R lidR package grid_terrain() giving weird plot

texts:
I am working on the pre and post lidar earthquake data from Kumamoto, Japan. I would like to generate a DTM with the grid_terrain(pre_online, res = 1, algorithm = knnidw(k = 10L, p = 0.5)) function in R Studio. When I plot the data with plot_dtm3d() I get a weird looking one. Is someone familiar with this problem?
The data is available at https://portal.opentopography.org/lidarDataset?jobId=pc1620409032095.

The plot looks like this and it should be a 3D terrain model.
pre_online
#> class : LAS (v1.2 format 1)
#> memory : 633 Mb
#> extent : -12000, -9000.001, -21500, -19500 (xmin, xmax, ymin, ymax)
#> coord. ref. : JGD2000 / Japan Plane Rectangular CS II (with axis order normalized for visualization)
#> area : 6 km²
#> points : 20.74 million
#> points density : 3.46 points/m²


texts: How to find a location is inside a country?

texts:
Given a Country Code and a Lat Lng pair, I want to find whether it is inside that country or not? What is the way to do it? Is there an algorithmic way to do this? I can accept certain amount of error.

I am node.js to implement this. Currently I do API call for Google to find the Country of the point. https://developers.google.com/maps/documentation/javascript/examples/geocoding-reverse But that doesn't full fill my performance requirement.


texts: Plotting multiple addresses with distance to specific POI

texts:
I'm not sure if this is the correct place for this question but here goes.

I am getting ready to move to a different city. Right now I have about 80 houses to choose from. One of my parameters for choosing a house is its distance from my work and from the church I plan on attending. What I would like to do is submit a list of addresses and get the distance/time between the house and my two POIs. Does anyone have any suggestions? I am a database developer/programmer by trade so I am not opposed to coding something if I need to.


texts: labelling polygon by clicking on polygon in QGIS

texts:
I would like to check if it is possible to have some sort of dynamic labelling of the polygon shapefile by clicking on selected polygon within this shapefile. Is this possible on the QGIS gui?
For example, I will have a base map + my generated sets of polygon layer. I would like to easily label this polygon layer by just clicking on polygons on the layer on the GUI. All the clicked polygon are labelled the same way.


texts: Merging points (bus stops) in QGIS manually or programmatically?

texts:
What I have:
I have a large data-set containing links (bus routes( (>2000) and their respective nodes(bus-stops) (>6000). The data set I received from the local bus operator is very messy. It contains multiple variations of names and multiple (nearby) locations for the same stop. There is no clear trend of these variations, so I cannot run a program (say, to merge all points within 50 m radius etc. Example in orange and yellow circle in the picture).

What I want to do:
I want to use QGIS clean the data manually.

Is there a method/plugin that lets you select multiple stops and merge their locations either to one of their locations, or even a manually entered location?

I now have ArcGIS (10.2) also, so a solution in either QGIS or ArcGIS is fine. Here is the screenshot from ArcGIS.

The points highlighted in yellow the circle are actual locations of the bus-stops. I want to merge them into one data-point(consider it as a node), while retaining their attributes. (attributes contain information about which routes pass through the stop and what order is the bus stop in a particular route etc ).

So if these data points are like this:

"Intersection xyz1 (Besides Public School)", lat1, long1, route 23, direction up
"Intersection xyz2 (Besides Govt. Hospital)", lat2, long2, route 23, direction down
"Intersection xyz3 (Besides General Pharmacy)", lat3, long3, route 87, direction up
"Intersection xyz4 (Besides Open Area)", lat4, long4, route 87, direction down

They should be converted to:

"Intersection xyz", lat, long, route 23, direction up
"Intersection xyz", lat, long, route 23, direction down
"Intersection xyz", lat, long, route 87, direction up
"Intersection xyz", lat, long, route 87, direction down


texts:
import Map from "esri/Map";
import MapView from "esri/views/MapView";
import Graphic from "esri/Graphic";
import GraphicsLayer from "esri/layers/GraphicsLayer";
import webMercatorUtils from "esri/geometry/support/webMercatorUtils";

async function kewlClickHandler(event) {

    function initDrawLine();

        function initDrawLine(Map, MapView, Graphic, GraphicsLayer, webMercatorUtils) {

          const map = view.map;

          var graphicsLayer = new GraphicsLayer(); // here
          map.add(graphicsLayer);

          ...... etc

No matter how I arrange my function params and header module includes; I repeatedly get the following error with my initDrawLine function:

Uncaught (in promise) TypeError: GraphicsLayer is not a constructor

New to arcgis-js-api, any ideas?

Update: This is the context, no errors or issues within stand alone CodePen:
in arcgis-js 4.13

But errors in actual environment in arcgis 4.11 - confusing.


texts:
I am trying to write a script partially in ArcGIS ModelBuilder and adding an if statement on my own.

It uses the Feature Compare tool in arcgis and I want to use the result.getOutput in an If statement so that if it is true then run a process - in this case a copy features.

The big picture here is to compare two features and if there are differences in them then truncate and update.

I need help with the If statement or alternative. I don't think I'm getting result.getOutput value to pass to the if statement and then run the if.  Can anyone see where this might be incorrect?

# -*- coding: utf-8 -*-
# --------------------------------------------------------------------------    -
# compare_test1.py
# Created on: 2015-11-03 14:06:44.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Local variables:
TONGIS_PARKING = "Database Connections\\ESRI10PUB GIS_VIEWER WAREHOUSE.sde\\TONGIS.PARKING"
TONGIS_PARKING__2_ = "Database Connections\\ESRI10PUB GIS_VIEWER WAREHOUSE.sde\\TONGIS.PARKING"
test = ""
PARKING_CopyFeatures = "G:\\Projects 10\\Corporate Services\\Information Technology\\Default.gdb\\PARKING_CopyFeatures"

# Process: Feature Compare
result = arcpy.FeatureCompare_management(TONGIS_PARKING, TONGIS_PARKING__2_, "GLOBALID", "ALL", "", "0.00001 Unknown", "0", "0", "", "", "NO_CONTINUE_COMPARE", test)
changes = result.getOutput(1)
#print "Comparision Result is: " + changes
change = str(changes)
if change == "True":
    # Process: Copy Features
    arcpy.CopyFeatures_management(TONGIS_PARKING__2_, PARKING_CopyFeatures, "", "0", "0", "0")


texts: Obtain elevation from a request lat,long with a map .tif

texts:
I am trying to get the elevation between two points in a route of Spain using APIs with HTTP request. I am using OpenStreetMap to get the routes, but I don't know how extract elevation with other tool.
I found out that SMRT has elevations from NASA and I have a .tif with the map, but I don't know:

how can I obtain the elevation from a .tif with a lat,long request in R?

Now I have this:

library(raster)
library(sp)
library(rgdal)

r<-raster("30n030w_20101117_gmted_bln075.tif")

test.dat <- structure(list(latitude = c(43.421604,43.42076,43.423283,43.422585,43.422546,43.21036), longitude = c(-5.825145,-5.8239202,-5.8272586,-5.8240676,-5.8240047,-5.822057)), .Names = c("latitude",
                                                                                "longitude"), class = "data.frame", row.names = c(NA, 6L))
points <- cbind(test.dat$longitude,test.dat$latitude)
pts2=SpatialPoints(points,proj4string = CRS("+init=epsg:4326"))
pts3=spTransform(pts2,projection(r))
coordinates(pts3)
extract(r,pts3)

With this code I obtain this values:

164 161 158 161 161 511

Are them the elevation in meters?


texts: Using Tobler's Function Equation in QGIS (Creating Friction Surface)

texts:
I'm struggling to get an equation working in QGIS that I've seen posted for ArcGIS. This article has a full description. Long story short, the author modifies the original Tobler's Hiking Function for use in ArcGIS, but it won't work in QGIS Raster Calculator. I've spent a few hours working on this with no luck.
The original formula is as follows:
(your_res /1000)/(6 * Exp(-3.5 * Abs(Tan(("your_slope" * 3.14159)/180) + .05)))

"your_res" = resolution/pixel size of DEM (in this case, ASTER)
"your_slope" equals slope raster (calculated using Raster -> Analysis -> Slope).
Here's how I modified the equation for QGIS Raster Calculator:
(your_res / 1000) / (6 * 2.71828 ^ (-3.5 * abs ( Tan ( ( "your_slope" * 3.14159)/180) + .05)))

The resulting raster returns few values and basically flattens the image out. Only the most severe slope changes show any variation. I don't have ArcGIS at home to test, but that article is fairly well known to archaeologists and I suspect it works.


texts: Error whilst exporting 3D pdf from Global Mapper

texts:
I’m try to export a 3D PDF from Global Mapper and I’m get this error message:

Error exporting to 3D PDF file.
Could not create U3D file.
PDF3DOverlay.cpp - 645 Version: v18.2.0 (64-bit)
Build Time: May 24 2017 14:41:24
Thread: Main UI Thread Windows 10 Home (64-bit) Memory: 3,636,162,560 of 8,018,649,088 available, GDI Usage: 425 GDI (Peak 535), 173 User (Peak 328)

The data I’m using is a 69mb .obj file
Does anyone have experience with this?


texts: Adding commas to number values in arc2earth generated KML balloons?

texts:
I'm using Info Window Builder in Arc2Earth to generate some kml for a large dataset.

I'm having some trouble formatting numeric values (pulled from the attribute table) to include commas (group separators).  Basically I'm trying to output the numeric values as US dollar amounts (1234 -> $1,234.00).

I've found some success using the floatformat filter as follows:
    [!MyFieldName|floatformat:2]
to successfully place exactly 2 digits to the right of the decimal point.

It seems I actually want to use the stringformat filter, but I can't figure out the correct syntax.
I've tried permutations similar to the following:
    [!MyFieldName|stringformat:("N", en-US)]
but with no luck.

Anyone know how to do this?


texts:
I have about 50,000 linestrings in PostGIS. They represent routes along roads. But I want them to be (roughly) pavements (aka sidewalks) instead. So I need to move all the lines 3 metres to the left of the direction they are travelling.

I found this function:

http://postgis.net/docs/ST_OffsetCurve.html

And tried to implement it, however some of my lines are self crossing, and thus it doesn't work. I think that this is reason anyway. An error I get it:

"ERROR:  GEOSOffsetCurve: IllegalArgumentException: Cannot get offset of single-vertex line"

I presume that this is becase there are some straight lines aswell as curves. So I need to cope with this.

Also, when I try to store these, sometimes it says that the result of st_offsetcurve has created a multi-string - though I'm not sure why?

Does anyone know of how I could do this please?

Thanks

James


texts: Summing volumes through a pipeline system

texts:
I need to fill a volume attribute for a pipeline system.  Volume would come from a well at the end and I can easily start with a spatial join to get volume for all the line segments at the edge of the system.

At this point I need to find a way to start summing the volumes through the system so that I know how much volume runs through a segment. Say a well volume 1 and well volume 2 come together. The resulting line needs to have a 3 applied to its volume attribute.

I understand that the best way to do this will probably involve building a geometric network and I have explored this but haven't come up with a good way to automate the summation.  I'd thought about using some of the hydrology toolset but had never worked with it before and don't know where to start.

Finally solved this using a combination of Attribute Assistant and the Water Utility Network Reporting Toolbar with the Summary Flow Accumulation tool.


texts: Can layer names in Layered PDF exported from ArcGIS for Desktop be overridden?

texts: Does anyone know how to override the layer names output by ArcGIS, where the output layer is called Other, or .  These usually refer to things like the map border, labels and the background.


texts: Getting rid of small objects on raster in ArcMap?

texts:
What I want to do is to get rid of 'dark' objects which are on the green areas.

My input looks like this:

And I want my result to look like here:

Raster value of green is 1, the rest has 0 value.
As You can see - I have the result, but i don't know wich tool in ArcMap I should use for this?


texts: Altitude point layer to slope heatmap in QGIS

texts:
I have a vector point layer (grid of 100km x 50km, a point every 20meters = 12.500.000 points). It has an altitude attribute.

I would like to have the information about the slope. I guess I need to compute the gradient but I have not found anything looking for "Gradient of point layer".
The final goal is to select zones given criteria like "slope is more than 10%".

How to compute the slope layer from the altitude layer?
Is heatmap ok to visualize and select from this slope layer?

I am very new to QGIS, so feel free to propose an alternative approach to this problem.


texts: Why can't I select all polygons in my DXF layer?

texts:
I am importing dxf file into QGIS. The file contains polygons. For some reason, when I select them only some of them are selected. Why can't I select all off them.

Here is how it looks like: (yellow is the selected polygons)


texts: Calculating shortest path using two network layers in PyQGIS

texts:
Using PyQGIS I can calculate the shortest path from point to point. But I need to calculate the shortest path based on two network layers: A and B. So that the result should be the shortest path from point to point following the lines of either layer A or layer B.
Here is the code I am using for calculating the shortest path using one network layer:
params = {'DEFAULT_DIRECTION' : 2, 'DEFAULT_SPEED' : 50, 'DIRECTION_FIELD' : None, 'START_POINT' : srcPointStr,
          'END_POINT' : destinationPointStr, 'INPUT' : QgsProcessingFeatureSourceDefinition('duct_b9e0ae0f_53bf_4875_9f30_e6e5e43060dc', True),
          'OUTPUT' : 'TEMPORARY_OUTPUT', 'SPEED_FIELD' : None, 'STRATEGY' : 0, 'TOLERANCE' : 0, 'VALUE_BACKWARD' : '', 'VALUE_BOTH' : '',
          'VALUE_FORWARD' : ''}
try:
   result = processing.run('qgis:shortestpathpointtopoint', params)
except QgsProcessingException:
   self.showErrorMessage("Error occures during Shortest Path Calculation.", "Error")

I have tried to pass an array of layers IDs as the INPUT, but it throws an error:
"Unable to execute algorithm Could not load source layer for INPUT: invalid value"


texts:
I have multiple cross-sections (point data).
The adjacent cross-sections are very far apart. I want to copy a set of cross-sections at a fixed interval (say n-meters) along a line. At this moment I'm manually doing it using Editor toolbar.

Figure below shows the before and after image.

Cyan are the original cross-sections and the red ones are the copied (copy of adjacent cross-section)

Is there any way to copy a set of points, at a fixed interval along a line perpendicular to it, in ArcGIS ModelBuilder?


texts: Looking for Alternatives to Google gmaps utility widget for ArcGIS JSAPI API

texts:
I have been using https://code.google.com/p/gmaps-utility-gis/ as a solution to displaying Google imagery as a base layer on a ESRI JSAPI map.  With the advent of JSAPI 3.4 and full AMD support this solution is showing it's age.

Is this still the best/only* solution for showing Google imagery on a JSAPI map?

*ESRI indicates that the WebTiledLayer will show Google imagery but does not provide an example.


texts: Using `CAST` and `ST_UNION` in same query in PostgreSQL

texts:
I want to union the geometry of waterlines grouped by a table column in a new view.
My problem is that I apparently needes to specify the geometry type using CAST function. But I cannot figure how to use the CAST function in combination with ST_UNION.

CREATE OR REPLACE VIEW _06_spildevand_vandloeb.drikkevandsledning_dissolve as
 SELECT row_number() over () as id,
    ST_UNION(the_geom)::geometry(MULTILINESTRING,25832),
    descr
  FROM _06_spildevand_vandloeb.ledninger_vand
  GROUP BY descr

This results in, though, that the name of geometry column is called ST_UNION, so I must have got the syntax wrong when using CAST and ST_UNION.

What is the problem?


texts: Can I increase the number of processes in ArcGis Server 10.1?

texts:
I use ArcMap 10.1. I have a map with many Mosaic and Raster layers. They are stored in many gdb-s. I publish this map with option "Create tiles on demand".

But this tiles are created very slowly (may be 2 minutes!!!). How increase the number of processes in arcgis server, which created this cache?

P.S. I know that ArcMap use only one сore. But maybe ArcGis server can use more cores?

Update

Thanks. I have the same problem. I set 24 cores. But CPU usages between 10-14%


texts: Why are some layers saved with custom CRS while others use appropriate EPSG code?

texts:
I have a number of .shp's that I have converted from EPSG 4148 to CPNO27. Some have come out perfectly and others revert to EPSG 1000007. I have checked and the only difference between the two is the addition of the bold text in CPO27 Cape/27.

As far as I understand, to save a file in a re-projection, one right clicks >> save as >> insert name >> ensure the EPSG is appropriately selected. So why would some of my .shp's revert to the 100007 EPSG (automatically might I add!!!) and not in CPNO27?

I have even gone as far as testing each file by setting CRS of .shop individually and if it changes spatial location, then doing a re-export, always selecting the CRS.

EPSG's follow:

1000007

+proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=0 +y_0=0 +a=6378249.145 +b=6356514.966398753 +towgs84=-136,-108,-292,0,0,0,0 +units=m +no_defs

Cape/27

+proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=0 +y_0=0 +axis=enu +a=6378249.145 +b=6356514.966398753 +towgs84=-136,-108,-292,0,0,0,0 +units=m +no_defs


texts: How to add a graphic table element in Python without Production Mapping extension?

texts:
I'd like to add a table in the Map Layout that includes only items visible at the current extent of the dataset. I know that the Production Mapping Extension lets you insert a Graphic Table Element filtered by the dataset current extent (I suppose). But we haven't this extension.

Has anyone ever tried to do something like this in Python?


texts:
I am not a QGIS expert.

I have to create a kml file to be used as a boundary.

I started using qgis with two different shapefiles that feature
I separated the polygons of the two shapefiles.
I merged the polygons that interested me thus creating two new shapefiles.
I converted the two shapefiles into two new kml files.
Since I did not succeed in merging the two kml files with qgis, I use an online kml merger and it worked.
I reduced the size of the newly merged kml file with qgis simplifying its geometry and then I filled the holes that formed inside the polygon. This kml file is well shown in Google Earth and in qgis.
Since I am finding problems in using this new kml file in the site where it should be used as a boundary, I decided to run the geometry check tool
the geometry check tool showed that the whole kml file has an invalid geometry (it is red in its entirety).

How do I fix this problem?


texts: Stop script working after two specific values are selected

texts:
I have a script tool that runs after values from multivalue box are selected.

The code is okay and runs withou any mistakes.
But is it possible to make a kind of exception for selection of fields GRAPH and ADVANCED_GRAPH? Because if I run script with two of these boxes activated I get lots of unnecessary data in  TOC.

I tried something like this:

all_selected = arcpy.GetParameterAsText(0)
selected = all_selected.split(";")
arcpy.AddMessage("Working...")

for s in selected:
    if s == "GRAPH;ADVANCED_GRAPH":
    arcpy.AddMessage("ERROR")

and this

all_selected = arcpy.GetParameterAsText(0)
selected = all_selected.split(";")
arcpy.AddMessage("Working...")

for s in selected:
    if s == "GRAPH" and if s == "ADVANCED_GRAPH":
    arcpy.AddMessage("ERROR")

Both methods don't work.


texts: Print layout of multiple maps revert to one identical map after save as image/zoom

texts:
I have a map with 8 graduated classes. I have created a print layout with 8 separate maps - one for each of the classes. I used "lock layers" to ensure that each separate map displays just one of the classes. Everything looks fine, however when I try to export the image (or PDF), all of the maps change to reflect the one I changed most recently and I end up with 8 identical maps. The image file does not save the different maps. I have updated to the most recent version of QGIS (3.10.13) in Windows 10.

Edit: it seems just zooming in or out of the print layout also causes this problem


texts: Converting total station points to cartesian coordinates

texts: I have a set of points which have been downloaded from a total station instrument. These data points are already in excel sheets and the issue now is how to convert these points (easting, northing and elevation) to cartesian coordinates in order to use for modeling purposes.


texts: Geojson in postgis Geometry : Z dimension error

texts:
i had this php code :

$geom1='{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-32514.400906806848,4810093.013501738],[-32877.47679116143,4809500.626532529],[-34473.099756614494,4810026.131101989],[-33412.535989157666,4811124.913383588],[-32514.400906806848,4810093.013501738]]]},"properties":null}';
$geom_poly=explode('Polygon',$geom1);
$geom_def=$geom_poly[0]."MultiPolygon".$geom_poly[1];
$geom2=explode('{',$geom_def);
$geom3=explode('}',$geom2[2]);
$CRS='"crs":{"type":"name","properties":{"name":"EPSG:3857"}}';
$geometria1=$geom3[0];
$geom4='{'.$geometria1.'},'.$CRS;

$uso_edit = $_POST['uso'];
$afeccion_edit = $_POST['afeccion'];
$retablo_edit = $_POST['retablo'];
$proteccion_edit = $_POST['proteccion'];
$altura_edit = $_POST['altura'];
$query = "INSERT INTO prueba(geom,prueba) VALUES  ( ST_SetSRID(ST_GeomFromGeoJSON('".$geom4."'),3857) , 'jsadfjlasd') ";
pg_exec($query);

I had this error :

pg_exec(): Query failed: ERROR: Geometry has Z dimension but column does not in

My PostgreSQL version is 9.4: and my table code is :

CREATE TABLE prueba
(
   id_0 serial NOT NULL,
   geom geometry(MultiPolygon,3857),
   id integer,
   prueba character varying(80),
   CONSTRAINT prueba_pkey PRIMARY KEY (id_0)
 )
 WITH (
 OIDS=FALSE
 );
 ALTER TABLE prueba
 OWNER TO postgres;

I try change the dimension with :

ALTER TABLE urbana_prueba
ALTER COLUMN geom TYPE geometry(MultiPolygon, 3857)
USING ST_Force_3D(geom);

but dosen't work......


texts: How to edit manually file geodatabase's OBJECTID column in QGIS?

texts:
I'm using QGIS version 3.4.1 with a Windows 10 OS. I'm trying to edit an Esri FGDB folder. Haven't succeeded. So far so I'm using a workaround: exporting specific layer to a new FGDB - with 'Save Feature As...' - and then editing.

This lets me edit but creates a new OBJECTID column. This column is sequential and does NOT correspond with the original table's OBJECTID.

How can I retain the old layer's fields schema (i.e. only one OBJECTID column) with the same values as the original imported layer?

Things I've tried:

Field calculator -

Leaving a blank space inside the 'Save Feature As...' FID option -


texts: How to indicate nodata into gdal_calc formula

texts:
I have a recurrent problem with gdal_calc.

Say I have two rasters I need to sum together, in which the same nodatavalue is set (e.g. --NoDataValue = 0).

In order to account for the fact that nodata + data = nodata, my formula would ideally look like this:

--calc="where(A==nodatavalue,0,A)+where(B==nodatavalue,0,B)"

Unfortunately, the "real" nodatavalue (i.e. zero) is not recognized (so anything like --calc="where(A==0)... would NOT work), nor I find any other way to refer to it.

Any ideas that DO NOT include the prior use of gdal_edit -unsetnodata?


texts:
In QGIS, I want to discern between a number of building types, based on topology. I have data on buildings, parcels and addresses.

I'm having a hard time to discern between row/terrace, end-of-terrace and semi-detached housing. Assume, for the sake of simplicity, that all other dwelling types (detached houses, appartment buildings, non-residential, etc) have already been classified and do not complicate things further. My working definitions are currently:

row/terrace: "at least 80% of opposite sides overlap with perimeter of another building; other sides do not overlap with perimeter of other building"
end-of-terrace: "at least 80% of one side overlaps with perimeter of a row/terrace building; other sides do not overlap".
semi-detached: "at most one side shares overlaps for 80% or more with perimeter of another building; if a second side shares overlap with the perimeter of another building, this overlap can be at most 50% of that side"

The latter ensures that buildings which are semi-detached in terms of living areas but that are attached to other buildings by their garages would count as semi-detached, to fit the local context of this project where such a house would be considered semi-detached.

So far so good. However, I do not know how to implement these rules in QGIS. How should I go about this?

How do I calculate perimeter overlap?
How do I distinguish between opposite sides of a polygon?
How does this work for non-rectangular buildings?

Is there an easier way to achieve what I want that I'm missing?


texts: KML Screen Overlay - works on google Earth but not on ArcGIS

texts:
I've been trying various things to get this to work on ArcGIS Online for the past couple of hours, but have not been able to figure it out.

It works perfectly fine on Google Earth but ArcGIS Online only shows the folder "Screen Overlay Folder" but not the image inside it.

Is there any specific syntax to use for ArcGIS Online?


texts:
I have multiple raster layers (TIFF) with names following this pattern:

"euj_LANDSAT7_250m_20170101_v2"
"euj_LANDSAT7_250m_20170102_v2">
"euj_LANDSAT7_250m_20170103_v2">
...

and I would like to rename them to follow this pattern:

LANDSAT7_170101.tif"
LANDSAT7_170102.tif"
LANDSAT7_170103.tif"
...

I've tried to use this code (edited according to PolyGeo's answer):

import arcpy
from arcpy import env

env.workspace = "E:/Training/TryArcpy/test01_2017"

in_data = arcpy.ListRasters()

for raster in in_data:
    out_data = raster[4:13] + raster[20:26] + fileExtension
    arcpy.Rename_management(in_data, out_data)

The message I get is:

Runtime error
Traceback (most recent call last):
File "<string>", line 3, in <module>
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py",
line 4537, in Rename
raise e
RuntimeError: Object: Error in executing tool

Is this because I have a raster list as in_data? If yes, how do I loop through multiple rasters without using a list, or use a list in appropriate way? And if no, what could be the problem?


texts: Duplicate source information from esri.tasks.PrintTask

texts: We are trying out the PrintTask in ArcGIS Server 10.1 pre-release.  When an Esri base map is enabled, we are getting source information in two places in the resulting .pdf.  I can confirm that this was not the case with 10.1 beta, and we are using the same .mxd's for the service.  Is there any way to remove one of these source references, or change the placement?  I didn't see a corresponding parameter in the javascript api.


texts:
I've just tried to set up Deegree on the server with Tomcat7 servlet by copying the WAR file to the webapps folder as described in the documentation. So now I can use the interface on the port localhost:8080/deegree-webservices-3.3.14. It is also indicated on the interface that the working directory is /usr/share/tomcat7/.deegree.
However, when I try to set the pw or load a workspace, I got the error notification that this folder doen't exist. Seemingly, no folder has been set anywhere on the server.

Should I do anything else to install Deegree properly?


texts: UnicodeDecodeError while opening file named by non-English chars in PyQGIS 3.X

texts:
I tried to open files named partially in Chinese via PyQGIS, the error appears:

Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/encodings/ascii.py", line 26, in decode
    return codecs.ascii_decode(input, self.errors)[0]
UnicodeDecodeError: 'ascii' codec can't decode byte 0xe5 in position 73: ordinal not in range(128)

Here is my code:

root = QgsProject.instance().homePath()
source_file = "/buffer_prac/1A27_八德店.csv?delimiter={}&xField={}&yField={}".format(",", "Response_X", "Response_Y")
uri = "file:/" + root + source_file

A27_lyr = QgsVectorLayer(uri, "A27", "delimitedtext")
QgsProject.instance().addMapLayer(A27_lyr)

The file can be opened manually by Layer -> Add Layer -> Add Delimited Text Layer. The file can also be opened properly while changing the file name into English.

I have checked solutions suggesting to change via sys.setdefaultencoding('utf8'),
however I later discovered that it is not encouraged to do so and it seems not working in my case either.

How can I open the file without changing filename via PyQGIS?


texts: How to calculate the number of points with differing characteristic in a polygon in ArcMap?

texts:
I have a .csv (exported to .shp) database with about two thousand entries of points. One column in this database has certain values (e.g. a point can be 'ash tree', 'beech tree', 'elm tree', etc). I want to visualise their relative counts within polygons (provinces in this case) in a pie chart. I'd also like to adjust the size of the pie chart based on the total sum in that province, but I suppose this will be easy.

Using ArcGIS 10.1 for Desktop.


texts: Adding tool as layer in table of contents of ArcGIS Desktop?

texts:
I have a custom script tool in my toolbox in ArcGIS Desktop 10.3. I want to add it as a layer in table of contents so when the map is shared people can have access to it easily. for example as shown in figure:

When I drag a tool from catalog to contents, it gives error:

  Could not add specified data object to the map

Can anyone help?


texts: Which expression are working in the HTML frame in QGIS print composer

texts: I'm using the HTML frame in the print composer. I would like to automatically insert data like $scale and the name field of the Atlas layer. Till now I only got the date working [%day($now)||'.'||month($now)||'.'||Year($now)%]. Is there a list of which expression are supposed to be working inside the HTML frame?


texts: Printing geological legend using Generic Mapping Tools?


texts:
As I ran several tools of sDNA I noticed these lines in every error report:

ERROR: Unable to obtain XYTolerance for near miss connection check
ERROR:   Either add a spatial reference to the feature class,
ERROR:   or override XYTolerance in the script environment settings.

At first I thought I have found where the mistake lies and filled the XY tolerance as the recommended value, yet everything remained the same.
(The projected coordinate system of this file is WGS_1984_EASE_Grid_Global, also in meters, so the two should be compatible, I guess?


texts: Error saving object from python script using GetParameterAsText

texts:
I have a script that takes in a line, and creates a certain polygon based on this.

My code obtaining the input information is as follows:

inFC = arcpy.GetParameterAsText(0)
hasClearway = arcpy.GetParameterAsText(1)
hasHighDivergence = arcpy.GetParameterAsText(2)
outFC = arcpy.GetParameterAsText(3)
fcName = outFC.rpartition("\\")[2]
fcPath = outFC.rpartition("\\")[0]
outFC = arcpy.CreateFeatureclass_management(fcPath, fcName, "POLYGON", inFC, "DISABLED", "DISABLED", inFC)

The polygon is finally created from an array and inserted into outFC

runwaypolygon = arcpy.Polygon(polygon_array)
cursor = arcpy.da.InsertCursor(outFC, ['SHAPE@'])
cursor.insertRow([runwaypolygon])

However I am getting the error 20469, "An Error Occured trying to save the object named X", when trying to set the output FC, see below screenshot:

I have looked into this http://support.esri.com/de/knowledgebase/techarticles/detail/20469
but this has not shed any light on the issue. Two further points:

(1) I have successfully used this method with an object of exactly the same type (a polygon produced in the same way).

(2) The script works fine when I ask the user to specify output file location and name separately, as here:

inFC = arcpy.GetParameterAsText(0)
hasClearway = arcpy.GetParameterAsText(1)
hasHighDivergence = arcpy.GetParameterAsText(2)
outFC_location = arcpy.GetParameterAsText(3)
fcName = (arcpy.GetParameterAsText(4))
outFC = arcpy.CreateFeatureclass_management(outFC_location, fcName, "POLYGON")

What is the error in my code or my approach?


texts:
I need to know how to plot a velocity field in QGIS 3.6.
I have u, v and magnitude data in a CSV file.

I installed some plugins such as Crayfish, vector field calc, contour plugin, point sampling tool to calculate velocity fields. But still I couldn't do it.
I found some tutorials on how to do it in QGIS 2.**. According to the tutorial I need Vector Field Renderer plugin as well. However, in QGIS 3.6 that plugin is not available. How can I proceed this further.

Problem solved using symbol library and some edits on csv files. https://docs.qgis.org/testing/en/docs/user_manual/working_with_vector/style_library.html#the-vector-field-marker


texts: How do I generate a random raster (or probabilistic cellular automaton) in QGIS?

texts:
I'm trying to create something like a cellular automaton land cover change model using QGIS, in which cells randomly change land cover based on probabilities determined from various input rasters.  I've already have a map of probabilities of conversion, and I need to take the next step.

I found this question, Randomly altering raster map of habitat types?, which looks very similar to what I want to do.  The suggested approach is to use a Markov chain or probabilistic cellular automaton, and the commenter writes: "The transitions can be implemented using Map Algebra statements, which makes this method practicable in any raster-based GIS, even those without direct or quick access to cell-level data."

Can anybody tell me how to do probabilistic map algebra in QGIS or one of its extensions?  I don't see any probability tools in the Raster Calculator.

Failing that, is there a tool that will simply generate random rasters to be later combined using the Raster Calculator?


texts: Will the index.html swf need to be re-compiled in order to embed the flexviewer app in between a html header and footer wrapper?  Or it this perhaps handled in the HeaderControllerWidget.mxml?


texts: Given city and state, how to get the nearest zip code in US

texts:
We have a situation where user enters a city, state for which we need to find a corresponding zip code. In case there is no corresponding zip code, we are required to find a zip code that close to the city, state.

Can you please throw some ideas as to how to go about doing this using Google maps/geocode API.

EDIT: Is it possible for US that a city entered by user does not have a zip code?


texts:
I have to develop an extension in C# for ArcMap 10.2 and I need to show information from a datatable in an a table window.

I can show the table without problems, but I have a problem I can´t resolve.

If I show table with one record, for example, and then I want to show the information of another table, even if the new table has several records, the table window show only one.

This is the piece of code I use:

if (pTable != null) { IStandaloneTable pStandAloneTable = new StandaloneTableClass();
pStandAloneTable.Table = pTable;
pTWindow.StandaloneTable = pStandAloneTable;
pTWindow.Application = Connector.ConnectorExtension.Aplicacion;

if(!pTWindow.IsVisible) pTWindow.Show(true);

if (pTable != null)
{
 IStandaloneTable pStandAloneTable = new StandaloneTableClass();
 pStandAloneTable.Table = pTable;
 pTWindow.StandaloneTable = pStandAloneTable;
 pTWindow.Application = Connector.ConnectorExtension.Aplicacion;

 if(!pTWindow.IsVisible)
  pTWindow.Show(true);

 System.Diagnostics.Debug.Print("Table: " + pTWindow.StandaloneTable.Name);
 System.Diagnostics.Debug.Print("Records: + pTWindow.StandaloneTable.Table.RowCount(null).ToString());

pTWindow.StandaloneTable.Table.RowCount has the correct number of records in the table but the table window doesn´t show them


texts: Plot Categorical Leaflet Polygons from a csv in R

texts:
I want to plot each polygon with the color stated in the column 'Categories' in the dataframe 'new'

This is what the dataframe 'new' looks like:

palcat <- colorFactor(new$Categories)

## to delete the special characters in each polygon value (new$Polygon)
new$Polygon <- chartr('()[]','    ', new$Polygon)

#subset Polygon column
blocks <- new$Polygon

####          PLOT TO LEAFLET          ####
map <- leaflet()

counter <- 1

for(i in blocks){
  m <- matrix(as.numeric(strsplit(blocks,",")[[counter]]),ncol=2,byrow=TRUE)
  palcat <- colorFactor(viridis(20), new$Categories)
  map <- addPolygons(map, data = m, color = 'black', weight = 0.25,
                     smoothFactor = 0.5,opacity = 1.0, fillOpacity = 0.1,fillColor = ~palcat(Categories))

  counter <- counter + 1
}

# map %>%  addPolygons(new$Polygon, fillColor = colorFactor(new$Categories))
map
map %>% addProviderTiles(providers$CartoDB.Positron)

If you change the color directly in fillColor = "blue", it works

    fillColor = colorFactor(new$Categories))

inside the addPolygon part but since it's inside the loop that turns each Polygon Cell to a Matrix its not working. Any ideas on how to fix this?


texts: Lost Polygons after clipping Shapefile using QGIS clip fuction

texts:
I'm trying to clip the CORINE Land Cover 10 ha Shapefile (https://gdz.bkg.bund.de/index.php/default/open-data/corine-land-cover-10-ha-clc10.html) with a district shapefile to get my investigation area. Thereby, I extracted the three relevant districts of my investigation from the original district shapefile (https://gdz.bkg.bund.de/index.php/default/verwaltungsgebiete-1-250-000-ebenen-stand-01-01-vg250-ebenen-01-01.html) with the following code:

import geopandas as gpd

districts_germany = gpd.read_file('C:/Users/lucas/Desktop/Master Thesis/Python/Data_New/District Data/vg250_01-01.utm32s.shape.ebenen/vg250_01-01.utm32s.shape.ebenen/vg250_ebenen_0101/VG250_KRS.shp')
districts_fr_em_br = districts_germany.loc[[203, 204, 205], :]
districts_fr_em_br.to_file('C:/Users/lucas/Desktop/Master Thesis/Python/Data_New/District Data/vg250_01-01.utm32s.shape.ebenen/vg250_01-01.utm32s.shape.ebenen/vg250_ebenen_0101/districts_fr_em_br.shp')

Then I used the clip function in QGIS with the clc10 shapefile (green) as input layer and the created districts_fr_em_br shapefile (pink) as overlay layer. However, after clipping, some areas of the clipped clc10 shapefile (yellow) are missing. Both shapefiles have the same Coordinate Reference System.

It seems to be that polygons, which go beyond the district boundaries, are missing. Is there a way to cut these polygons at the boundary? However, I think also polygons, which are not crossing the border, are missing. Thus, I'm not sure if this is the problem.
I also tried to convert the MultiPolygons to Polygons and clip it afterward. But the same areas are missing.


texts: LiDAR Basics - How classes are determined at LiDAR data?

texts: LiDAR data have classes as Medium vegetation, Buildings, Ground, Model Key and more. What determines these classes ? is it combination of intensity and scan angle rank ? how to reclassify these classes using ArcGIS ?


texts: Converting raster into geostatistical layer?

texts:
I have downloaded a raster surface (SRTM) with pixel size equals to 30*30, and I want to use that surface as a geostatistical layer for further analysis. Is there any available tool that converts directly raster surface into geostatistical layer? I am using ArcGIS 10.4.1.

Therefore, if it does not exist any available tool, this methodology sounds accurate?
1) Considering the extent of the raster surface, I will construct grid points every 30 meters.
2) Using extract values to points tool, I will attach Z coordinates as attributes of grid points layer.
3) Interpolation of grid points using IDW tool(geostatistical analyst toolbox), gaining a new updated SRTM.

How accurate is this procedure comparing the raw SRTM with the modified SRTM that was developed through grid points?


texts:
I'm working on a student project trying to make a slope map with degree values. I downloaded a DEM file from a map service and used the raster ->analyze(tool)->terrain model and got the model. Now I would like to make a slope map

I googled it and found Calculating slope in percentage with QGIS

The second answer seem to have the solution, but I'm simply too much of a beginner to quite get what I have to click, step by step, to do this. He links to gdaldem

How do I make the operation to get the percent slope values?


texts: Pathfinder Office - 0% base data coverage

texts:
I am trying to post-process data collected with a Trimble Juno SB (sometimes hooked up to a ProXRT), in Pathfinder Office version 4.20.  I had no issues when I tested this in the US.  But now that I am in the field (Iceland), I keep getting " 0% total coverage.  No matching base data found."

The same files are working just fine on my colleague's computer.  Therefore, there must be something wrong with my configuration, but I can't figure it out - there are no obvious differences between my software setup and hers.  Any suggestions very welcome.

Here is the correction log.  Note that the Local Time for the base data spans 0 seconds, even though the reyk1740.15d file clearly contains data for the entire day.  My log files from data collected in the US, UK, and Ireland include a 24 hour span on this line.  What's wrong with my base data configuration, specific to my computer with data collected in Iceland??  Any suggestions very welcome!

I have the same problem regardless of which Icelandic base station I choose.

Searching for base files...

File [path redacted]\Iceland2\Base\SOPAC, Reykjavik, daily\reyk1740.15d.Z downloaded.
File [path redacted]\Iceland2\Base\SOPAC, Reykjavik, daily\auto1740.15n.Z downloaded.
Unable to transfer files.
Unable to download file from url ftp://garner.ucsd.edu/pub/nav/2015/174/auto1740.15n.g.
Successfully found or downloaded 2 of 3 files.
Search complete.

--------Base Data Details:--------------------

Using reference position from base provider:

Name: SOPAC, Reykjavik, daily
  Position: 64°08'19.61761"N,  21°57'19.74233"W,  93.05 m
Source: [path redacted]\Iceland2\Base\SOPAC, Reykjavik, daily
reyk1740.15d.Z
    Local time: 6/22/2015 23:59:46 to 6/22/2015 23:59:46
    Position: 64°08'19.61688"N, 21°57'19.74045"W, 93.14 m, 0.06 m Antenna height
    Distance from base provider: 0.03m

--------Coverage Details:--------------------
Rover file: [file1]
  Local time: 6/23/2015 12:58:35 to 6/23/2015 14:39:46
  0% total coverage.  No matching base data found.

Rover file: [file2]
  Local time: 6/23/2015 14:40:30 to 6/23/2015 16:31:50
  0% total coverage.  No matching base data found.

Click Confirm to continue...


texts:
I have an input csv file that I process using a python script and generate a new output csv file. I then, in the same python script use Table to Table (Conversion) in order to store the data from the output csv file in a file geo database.

The information in the input csv file will change but I'm currently only able to create a new and not update the table. So the question is, how can I update the table or if there's a smarter solution to this problem than generating and converting csv files?


texts: Euclidean Distance in Model Builder error

texts:
I'm trying to use Model Builder to create Euclidean Distance raster files for a set of line features.

Each line also have a buffer polygon that corresponds to each line that I want to have as an extent. However, it doesn't seem to be a way to iterate extents and I get a different looking result file when I use the model in comparison to using only the tool.

How can I run the model so it generates one unique Euclidean raster file for each and everyone of my lines while having the correct buffer as extent? I can do it manually but that would take forever for all of the lines.

Edit:

I made some progress in getting some results. However, I would need some way to iterate the extent each time. If I put in no extent the model crashes and if I use one buffer polygon for extent everything is in the same place, which isn't what I want to achieve.


texts: Separating Overlapping DMAs by State using PostGIS?

texts:
I have two shapefiles. One is the DMA and another of states

I know that some DMAs, such as Chicago will overlap between multiple states, for instance the Denver DMA overlaps states outside of CO and into UT.

Is there a way to recreate DMAs by states, so it would be something like Denver DMA - CO and Denver DMA - UT (to differeniate which states the DMA belong to)? Basically I'm thinking of a way to draw the state boundaries into the DMA maps and separating it accordingly.


texts:
I have obtained a desired value from a dataset using:
ogrinfo input.shp -dialect SQLite -sql "SELECT SUM(field) from input

The ogrinfo command prints the desired value via the command line when I run this in a python script. I need to write out this result to a .txt or .csv file as my intention is to display it in a barplot, ideally using plotly.

Other relevant info:
Windows 7, Python 3.3
GDAL 1.11


texts: Is there an algorithm to perform a smoothing of a linestring and maintain the original points?


texts:
I am trying to read a hyperspectral image in Uint16 .tif format and visualize it.
I use GDAL in Python to read the data and everything works perfectly when I use matplotlib to show one band of the data (The first image). The image however becomes completely corrupted when I want to concatenate(stack) 3 different selected bands and visualize them. Here is the code:

import numpy as np
import struct
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo.gdalconst import *
import matplotlib.pyplot as plt

cube = gdal.Open(filename)
bnd1 = cube.GetRasterBand(29)
bnd2 = cube.GetRasterBand(19)
bnd3 = cube.GetRasterBand(9)

Now I turn each band into a ndarray:

img1 = bnd1.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)
img2 = bnd2.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)
img3 = bnd3.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)

and stack them to have a 3 band image

img = np.dstack((img1,img2,img3))

f = plt.figure()
plt.imshow(img)
plt.show()

here are the results:


texts:
I am new to Python and I am using the Python console in QGIS 3.10 to make the coordinate grid on a vector layer. I had lots of vector files which are basically square shaped of 500 meters but all are of different location. I want to export this vector file but after making grid on top of that. I want to do it in Python as I can loop it and do all the process at once. My workflows goes like this:

Load vector layer
Open map canvas
Display vector layer in map canvas
Make coordinate grid on top of that
Export that as image

Here is my code and in output it is nit creating the grids
am new to Python and I am using the Python console in QGIS 3.10 to make the coordinate grid on a vector layer. I had lots of vector files which are basically square shaped of 500 meters but all are of different location. I want to export this vector file but after making grid on top of that. I want to do it in Python as I can loop it and do all the process at once. My workflows goes like this:
Load vector layer
Open map canvas
Display vector layer in map canvas
Make coordinate grid on top of that
Export that as image
Here is my code and in output it is nit creating the grids
fn = r"C:\Users\Asus\OneDrive\Documents\map secien\360m\parent
> grid\grid_no_A1.shp" layer = iface.addVectorLayer(fn, '', 'ogr')
>
> from qgis.PyQt import QtGui   layers =
> QgsProject.instance().mapLayersByName('grid_no_A1') layer = layers[0]
> project = QgsProject.instance() manager = project.layoutManager()
> layoutName = 'Layout1' layouts_list = manager.printLayouts()
> # remove any duplicate layouts for layout in layouts_list:
>     if layout.name() == layoutName:
>         manager.removeLayout(layout) layout = QgsPrintLayout(project) layout.initializeDefaults() layout.setName(layoutName)
> manager.addLayout(layout)
> # create map item in the layout map = QgsLayoutItemMap(layout) map.setRect(20, 20, 20, 20)
>
> # set the map extent ms = QgsMapSettings() ms.setLayers([layer]) # set layers to be mapped rect = QgsRectangle(ms.fullExtent())
> rect.scale(1.0) ms.setExtent(rect) map.setExtent(rect)
> map.setBackgroundColor(QColor(255, 255, 255, 0))
> layout.addLayoutItem(map)   map.attemptMove(QgsLayoutPoint(5, 20,
> QgsUnitTypes.LayoutMillimeters)) map.attemptResize(QgsLayoutSize(180,
> 180, QgsUnitTypes.LayoutMillimeters))
>
>
> #grid = QgsLayoutItemMapGrid(layout) out = QgsLayoutItemMap(layout) grid = QgsLayoutItemMapGrid("hr", out ) grid.annotationEnabled()
> grid.setIntervalX(0.1) grid.setIntervalY(0.1)
> grid.setAnnotationEnabled(True)
>
>
>
> layout.addLayoutItem(out)
>
> layout = manager.layoutByName(layoutName) exporter =
> QgsLayoutExporter(layout)   fn =
> r'C:\Users\Asus\OneDrive\Documents\map secien\layout_export.jpeg'
> exporter.exportToImage(fn, QgsLayoutExporter.ImageExportSettings())


texts: How to disable "Simplify geometry" via Python?

texts:
After I've loaded a layer from a table in my db I'd like to disable the 'Simplify geometry' (found in Layers/Properties/Rendering) programatically using Python.
I've tried several versions of QgsVectorSimplifyMethod.SimplifyHint

current_layer = iface.mapCanvas().currentLayer()
if current_layer:
  try:
    current_layer.setSimplifyMethod(QgsVectorSimplifyMethod.setSimplifyHints(QgsVectorSimplifyMethod.SimplifyHint.NoSimplification))
  catch Exception, e:
      print e

but can't get the usage of QgsVectorSimplifyMethod.SimplifyHint.NoSimplification correct.


texts: Are there any NTS index maps showing the most recent update?

texts:
For the Canadian National Topographic System, are there any index maps available that do not only show map numbers throughout the country, but also the most recently available update for each map?

As a bonus, is there any publicly available list showing the plan for updates for the next several years, for NTS maps?

I could not find either on the NTS website.


texts:
I have spatial data that has been put together with ArcGIS for Desktop and exported to me as a .pmf file.

How can I use the .pmf file in QGIS?

There are about 1000 points and I don't want to have to redo this work.

Do I have to convert it to something else so I can import it to QGIS? If so, can someone suggest a converter program?

I also have the points in an Excel file.


texts: Setting color of feature depending on attributes with PyQGIS?

texts:
I'm trying to set the color of a vectorlayer feature depending on the feature attributes with PyQGIS.

How do you acces the attributes of a vectorlayer feature?