texts: KML Layers suddently not appearing on embedded Google map

texts:
I have an embedded map on:

http://wtc.army.mil/modules/family%20and%20caregivers/sfacMap.html

I am receiving the error below in different forms related to the session id/token trailing after the GetViewportInfo? part of the request whenever I zoom in/out.

NetworkError: 500 Internal Server Error - http://maps.googleapis.com/maps/api/js/ViewportInfoService.GetViewportInfo?1m6&1m2&1d-52.02177634131835&2d-180&2m2&1d90&2d180&2u3&4sen-US&5e4&6sr%40275000000&7b0&8e90&9b0&10b1&callback=xdc._s8nn5v&token=34036

This all of a sudden occurred when previously my marker overlay was working. The functionality was that the markers would display info from my kml when clicked. Currently this does not seem to be working and I have made no changes since it was functional.


texts: Mapping origins and destinations to the generated isochrones in R

texts: I created isochrones in QGIS using Iso-area as polygon algorithm with an interval of 30 minutes. I exported isochrones as a geopackage file and read it into R using st_read() function. I also have some coordinates related to some origins and destinations. I need to map these origins and destinations to the created isochrones to identify what percentages of them lie within that 30 minutes interval. In other words, if I want to travel from the origins to the destinations what percentage of those trips lie completely within the created isochrones? Is there any package in R to help with achieving this goal?


texts: Is there QGIS Portal similar to Portal for ArcGIS?

texts:
I'm a beginner.

ArcGIS Pro provides a wide range of maps/data via its "All Portal". They have a few excellent maps I would like to use.

Is there an equivalent service from within QGIS or do I need to subscribe to ArcGIS to utilise these maps?


texts: How to get ESRI Addin to read the System.ServiceModel configuration?

texts:
I have a complex set of add-in's that we need to run in ArcMap, ArcServer, and ArcGIS Explorer for Desktop.  We are experimenting with changing from the old Web Reference service connection to a WCF Service Reference connection.

When I try to create the client it cannot find the configuration settings to create the client.  When I add just the URI to the embedded Settings.settings file, and create an EndpointAddress and Binding object to pass to the client, I get a client but my first attempt to access the server fails.

So how do I access the WCF service configuration data from inside an ArcMap Addin?


texts:
I am trying to read NetCDF4 files from data developed by Toshicika Iizumi (https://doi.pangaea.de/10.1594/PANGAEA.909132) on global crop yields. I would like to turn the data into a form that is readable by ArcGIS Pro since ArcGIS doesn't recognize the NetCDF4 format.
I am reading the file with xarray like this...
import rioxarray
import xarray
fn = ".../yield_2012.nc4"
xds = xarray.open_dataset(fn)
xds

There doesn't seem to be any information on the projection.
How can I turn this into a raster or GeoTIFF for use in ArcGIS Pro?


texts: Passing field variable as distance with Drive Time Trade Areas using ArcPy with ArcGIS Pro

texts:
I have a dataset of 600,000 records where I need to calculate a drive time that is different for each based on model that estimates the actual/variable drive time needed at that location.  Within ArcGIS Pro (2.5) I can create a drive time trade area but it only allows input for 1 specific drive time point.

How can I change the python code so that it reads from the record/file the specific drive time that I need for it to calculate?

Here is the Business Analyst drive time code:

arcpy.ba.GenerateDriveTimeTradeArea("Heat_Map_Grids_Point", r"C:\Default.gdb\Heat_Map_Drive_Times", "Driving Time", [5], "MINUTES", "GRID_ID", "OVERLAP", "KEEP_OVERLAP", "TOWARD_STORES", None, "TIME_ZONE_AT_LOCATION", None, "STANDARD")

the [5] (the distance value parameter) is the variable I need to convert to a field value.

Is this even possible, and what would some example code look like?

I am thinking sending the field values to a list and looping through that list would do the trick...but not sure how to create that code.


texts: Seeking vector data of agricultural field boundaries in small area of south west England?

texts:
I would like to find a vector polygon dataset that shows field boundaries, along with tracks, roads, paths etc, for a small area of south west England - about 5km by 5km.

Are there any free datasets available?

The OpenStreetMap data is from the 1940s and out of date.


texts: What free software can be used for print-quality map production?

texts:
I am creating art prints using spatial data.

I use ESRI software every day for my regular job, therefore I am proficient in GIS. However, I can't afford ESRI software for my art projects at home.

Currently I am using QGIS. What are some recommended products for producing print-quality maps?


texts: How to add a script tool in a toolbox without the Add Script Wizard?

texts:
I have been put on a project where I need to modify and optimize an existing PythonToolbox and associated Toolbar (each button calls a script). I want to add a new script tool in there but the Add Script Wizard isn't available on my Basic License (should it? when I right click on the toolbox and go to Add, the Script option is disabled).

What I've done is add a new class in the .pyt file and updated the self.tools of the Toolbox initialisation.

I've also, and I'm pretty this isn't the right thing to do, created a .xml file for the new tool with Notepad++. Basicaly I copied the file from another tool in the toolbox and modified the file name as well as the tool name under <resTitle> and <tool name="" displayname "">.

And of course what I get is a tool that is recognized but invalid.

What would be the correct procedure to add a script tool outside of the Add Script Wizard? How should I handle the new .xml file so it links with the .pyt file?


texts: Is there a way to publish map server by command on ArcGIS Server 10.1

texts:
I create a service define file by ArcMap10.1, but it's very large(about 90G).

Now I want to publish it to Server by Server Manager(access Brower).

Unfortunately, the manager is always time out (maybe the sd file is very large).

Is there another way to publish the map server?


texts: Error code while downloading cloud free mosaic Sentinel-2 image in QGIS

texts:
I'm trying to create a composite image of 2 Sentinel-2 images to make it cloud free using the semi-automatic classification plugin on QGIS, but when I press 'run' in 'mosaic band sets' in the plugin I receive the following error message (also I am new to QGIS and have never used Python before).
How can I resolve this issue?
2020-12-21T19:48:47     WARNING    Traceback (most recent call last):
  File "C:/Users/heidi/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\maininterface\mosaicbandsets.py", line 69, in mosaicAction
  self.mosaicBandSets()
  File "C:/Users/heidi/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\maininterface\mosaicbandsets.py", line 228, in mosaicBandSets
  File "C:/Users/heidi/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\core\utils.py", line 6912, in addRasterLayer
  for i in range(len(mcList)):
 AttributeError: 'NoneType' object has no attribute 'setName'


texts: Splitting lines at certain lengths in ArcGIS ModelBuilder

texts: How to split lines at a certain length in the ArcGIS ModelBuilder using below the split tool from the Editor Toolbox?


texts:
I am trying to strip data out of an existing shapefile to create new shapefiles. I have made the loop work, and made the new file, but I cannot work out why I cannot add records to the new file I am writing:

import shapefile as shp

coord1 = (39.737094, -121.543857)
coord2 = (39.651237, -121.403610)

sfinput = shp.Reader(r'C:\...\fire_archive_M6_98806.shp')
w = shp.Writer(r'C:\...\Output')
w.fields = sfinput.fields

for i in range(len(sfinput.records())):
    #test if within latitude
    if sfinput.record(i)[0] >= coord2[0] and sfinput.record(i)[0] <= coord1[0]:
        #test if within longitude
        if sfinput.record(i)[1] <= coord2[1] and sfinput.record(i)[1] >= coord1[1]:
            w.record(sfinput.record(i))

The error I am recieving:

record = [recordList[i] for i in range(fieldCount)]
IndexError: tuple index out of range


texts:
In PyQGIS can I get the total length of a road where it's broken up into segments?

In SQL I would run

SELECT
 format_number(ROAD_ID,0) as Road_ID,
 ROAD_NAME as Road,
 format_number(Sum(st_length(geometry)),2) as Total,
 Count(*) as Segment_Count
FROM
 CSC_Roads
GROUP BY Road, Owner
ORDER BY Owner, Road

For each segment I do have GIS_Length which is just $length but how would I do a aggregate sum when the RoadMntnc name is the same so that for each road segment we have the Total_length and Segment_length.

Would the best solution be something like?

geom=f.geometry()
ntext=to_string(format_number(sum(geom.length(),group_by:="RoadMntnc"),0))

I keep getting invalid syntax on the above

In field calc the following works - to_string(format_number(sum($length,group_by:="RoadMntnc"),0))

==== The complete code is

import re

def Update_Field(fieldName, ntext): # Function to simplify updating fields
    print (ntext)
    fieldIndex = layer.fields().indexFromName(fieldName)
    dpr.changeAttributeValues({f.id(): {fieldIndex: ntext}})
    layer.commitChanges()

RoadType_dict = {'Rd':'Road', 'Avenue': 'Ave'}
layer = iface.activeLayer()
selection = layer.getFeatures() #all features
#selection = layer.selectedFeatures() #only selected features

#Locality Data
Locality_layer=QgsVectorLayer(r'Z:\My Drive\Mangoesmapping\Spatial Projects\Gen_Library\DSC\vector\DSC_Localities_areas.shp', "Locality", "ogr")

dpr = layer.dataProvider()

for f in selection:
    print (f.id())
    text=f['RoadMntnc']
    print (text)
    stext=re.split('_+', text)
##Update RoadName
#    ntext=' '.join(stext[:-1])
#    fieldName='RoadName'
#    Update_Field(fieldName,ntext)
##Update Road Type Abbreviation
#    ntext=stext[-1:][0]
#    fieldName='Type_Abb'
#    Update_Field(fieldName,ntext)
##Update Road Type
#    ntext=RoadType_dict.get(stext[-1:][0])
#    fieldName='Type'
#    Update_Field(fieldName,ntext)
##Update Road Length
    geom=f.geometry()
    ntext=to_string(format_number(sum(geom.length(),group_by:="RoadMntnc"),0))
    fieldName='Chainage'
    Update_Field(fieldName,ntext)

layer.commitChanges()
layer.updateFields()


texts: lang.hitch use in ArcGIS API for JavaScript

texts: Can someone explain lang.hitch() function use in ArcGIS API for JavaScript?


texts: Copy files to a new folder if the file name matches value from a specific field in a feature class

texts:
I have this folder with bunch of .jp2 files, and from this folder I need to copy only the files that their name matches the value name from a specific field in a feature class. Right now I can print the values name, but kind of stuck on the copy function. I don't get an error when running the script, but nothing gets copied.

What am I doing wrong?

import os, shutil, arcpy

workdir = os.path.join(r'G:\County','JPEG2000','1Foot')
dst = os.path.join(r'T:\gis','Images2017','jp2','1foot')

allFiles = os.listdir(workdir)

path = r"T:\gis\airphoto\Images\Images2017_Mosaic.gdb"
arcpy.env.workspace = path

fc = r"T:\gis\airphoto\Images2017\Images2017_Mosaic.gdb\Extent_1ft"

rows = arcpy.da.SearchCursor(fc)

fields = arcpy.ListFields(fc,'FILENAME')

for row in rows:
   for field in fields:
      fieldName = field.name
      value = row.getValue(fieldName)
      print (field.name+": "+ str(value))

for aFile in allFiles:
   if aFile == str(value):
      shutil.copyfile(workdir,dst)


texts: Displaying inner area and removing outside area when clicking on layer in OpenLayers 3?

texts:
I want to display multiple WFS layers on OSM layer.

If i selected any one layer from multiple layers then display only that layer inside area and cut or remove the outside area of selected layer.

Can anyone give me sample example to do this in OpenLayers 3?


texts: Layer stacking Multispectral reflectance imagery in ERDAS Imagine 2016

texts:
I have Green, red, Red edge and Near infra-red reflectance maps captured with a Parrot Sequoia camera and post processed with Pix4D.
They are all Single Float Tiff files, around 500Mb in size.

I would like to Layer stack them using ERDAS Imagine 2016.

The individual band images display correctly (image 1), but when i stack them using the spectral Layer stack tool, the resulting image does not seem to display the reflectance data correctly, with each band layer having one block colour (image 2).

Applying a stretch to this image can pick out more information, but it's far from the original imagery.

Image 3 shows the method i have been using to stack the imagery. I select band 1 from each individual image and add them in Green, Red, Red-edge, NIR order.

Does anyone have any idea what i could be doing wrong?


texts: Automating capture of Extent from multiple mxds with multiple dataframes

texts:
I want to create a single shapefile from multiple mxd's that have multiple frame sets with different extents in them. I have found/started a script to do this (attached) but can't figure out how to write the captured X&Y Max/Min into the shapefile that is created for this. See output below. I also want it to write the scale and title of the frame as well as the file name of the mxd.

Can you help complete this script?

 import arcpy, glob, os from arcpy
 import env from arcpy import mapping
 env.overwriteOutput = True

 path = r"C:\temp" mxdList =
 glob.glob(path + "\*.mxd")

 env.workspace = r"C:\temp\Test.gdb"

 y = 1

 for mxd in mxdList:
     mxd2 = mapping.MapDocument(mxd)
     dataframe = mapping.ListDataFrames(mxd2, "*")[0]
     frameExtent = dataframe.extent
     XMAX = frameExtent.XMax
     XMIN = frameExtent.XMin
     YMAX = frameExtent.YMax
     YMIN = frameExtent.YMin
     pnt1 = arcpy.Point(XMIN, YMIN)
     pnt2 = arcpy.Point(XMIN, YMAX)
     pnt3 = arcpy.Point(XMAX, YMAX)
     pnt4 = arcpy.Point(XMAX, YMIN)
     array = arcpy.Array()
     array.add(pnt1)
     array.add(pnt2)
     array.add(pnt3)
     array.add(pnt4)
     array.add(pnt1)
     polygon = arcpy.Polygon(array)
     arcpy.CopyFeatures_management(polygon,
 "Polygon_Extent" + "_" + str(y))
     y = y + 1

 list = []

 lstFCs =
 arcpy.ListFeatureClasses("Polygon_Extent*")
 for fc in lstFCs:
     list.append(fc)

 arcpy.Merge_management(list, "Extent")

 for item in list:
     arcpy.Delete_management(item)

Even with this I get errors...

Traceback (most recent call last):
  File "P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\extent_creation.py", line 32, in <module>
    arcpy.CopyFeatures_management(polygon, "Polygon_Extent" + "_" + str(y))
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 1943, in CopyFeatures
    raise e
ExecuteError: ERROR 000210: Cannot create output Polygon_Extent_1
Failed to execute (CopyFeatures).


texts: Calculating area weighted average in ArcGIS Desktop?

texts:
I want to calculate area weighted average of districts of India.

I have Particulate matter 2.5(PM2.5) data for Districts of India, and now I want to find area weighted average of PM2.5 for each of these districts.


texts: How to find the ID of the selected object in the layer? PyQGIS 3.0

texts:
How can I get a list of the ID of these selected objects?

QGIS 3.0


texts:
In QGIS Python, how to get alias of field. In class QgsField, i can not find method to do this.

lay = iface.activeLayer()
provider = lay.dataProvider()
fields = provider.getFields()

for field in fields:
    print field.[SOMETHING]


texts: Iterating over line features to change colors

texts:
I want to iterate over line features to change color of each feature to the color I choose.
So my file has 7 features, I want am writing the code below, but it does not seem to work.
line_feats = [f for f in line_layer.getFeatures()]
clr = ["blue","red","green","yellow","magenta","cyan","black","red"]
for l in line_feats:
    line_layer.select(l.id())
    iface.mapCanvas().setSelectionColor( QColor(clr[l.id()]))
line_layer.triggerRepaint()

I have tried printing l.id() and it goes from0 to 7, so indexing is not the issue. The thing is, it does not change the color individually, rather it changes for the entire file.


texts:
I am using AcadMap to bring in .shp files for polylines. I have AutoCAD tickmark block at midpoint but it is at 0 degrees (straight up).

Is there any way to rotate block angle to match polyline angle?

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Revised 1:52 pm

The .shp file contains lines representing electric power lines which has object data attached. (Length of line, material, etc.)

The block is in AutoCAD with attributes. In the MAPIMPORT dialog box under Point Mapping you can tell esri to use certain blocks to represent esri symbols. The block comes in ok, but it is not aligned with line.
(See attachment).

I would like for esri to rotate block to align with line.


texts: How can I round to certain # of decimal places in data-driven field display?

texts:
Can I round the value being shows in a textbox that is provided through Data Driven Pages?

In the example below I have a SHape_Length value I want to display on the map but even though I have sent the number of decimal places to show in the table, that is not being transferred to the table.


texts: Adding Direction and Distance into Attribute table

texts:
I just want to add direction (Bearing: i.e N 25 35 E) and Distance (Length: 125 meters) as my new field in polyline/line data.
Is there a Plugin to generate these fields?
I tried to used "Export/Add geometry columns" in my line data, but only "Length" value was added.


texts: QGIS plugin + project: Enter db credentials(password and username) programmatically

texts:
I am developing QGIS Python plugin. I have a QGIS project and layers stored in a PostGIS database.
I don't want to store the password into the project and I want to avoid that user need to type credentials, and also I want to avoid storing the credentials permanently.

When I open the project, I'm immediately prompted for credentials before any layers from the db are rendered.

Is there any way to enter db credentials(password and username) or set datasource programmatically before dialog is shown?

I want avoid this dialog to be shown:


texts:
So I've got a number of different rasters with different values that represent similar data. I'm trying to create a default styling for them in geoserver. I've imported the QGIS style into geoserver. However the issue arises where the value ramp is based on the original raster's values.

How do I change the SLD so that it may be used for multiple rasters and not just the original?

Relevant SLD Code:

                    <sld:ColorMap>
                        <sld:ColorMapEntry color="#000000" label="Min" opacity="1.0" quantity="0"/>
                        <sld:ColorMapEntry color="#ff0000" label="Half" opacity="1.0" quantity="126"/>
                        <sld:ColorMapEntry color="#ffffff" label="Max" opacity="1.0" quantity="255"/>
                    </sld:ColorMap>


texts: Output of layer consisting pixel value 1 and 0 totally black?

texts:
I am computing the NDWI of an area , and I used the code below to generate a result of pixel value 1 and 0 only. Pixel value 1 will be displayed in blue means the pixel's NDWI value once had exceeded 0.3 and pixel value 0 will be displayed in white meaning the pixel's NDWI value has never reach 0.3 or more

var gte3 = ndwi.map(function(image){
  return image.gte(0.3);
});
var ndwigt3 = gte3.max();
var ndwigte3Viz = {min: 0, max: 1, palette: ['white', 'blue']};
Map.addLayer(ndwiget3, ndwigte3Viz , 'ndwigte3', false);

I could see the white and blue pixels when i open the layer in code editor, but once i output it with the following code, it turns out to be a total black image

Export.image.toDrive({
  image: ndwigte3,
  description: 'ndwigte3',
  region: Border,
  scale: 30,
  maxPixels:1e13,
  skipEmptyTiles: false
  });


texts: Importing existing data model into geodatabase

texts: How can I integrate an existing ArcGIS datamodel like Arc Marine into my geodatabase? I don't have schematics in my ArcGIS 9.3.1.


texts: Can you rename symbol once it’s saved in style Reference in ArcMap?

texts: Is there a way to rename a symbol in Style Reference file without deleting the file and starting all over again?


texts: Getting difference between field and constant value in ArcPy?

texts:
I need to calculate the difference between two input parameters within a table. the first parameter is a field in the table, and the second is a constant i want to subtract from this field, both of which are asked for in a script/tool.

currently my code looks like this:

heightfield = arcpy.GetParameterAsText(4) # field (float) in table inObsPtswHeight
heightOffset = arcpy.GetParameterAsText(5) # linear unit, constant value
arcpy.CalculateField_management(inObsPtswHeight, heightfield, heightfield - heightOffset, "PYTHON")

This is run in an arcgis script toolwindow and returns the error:

arcpy.CalculateField_management(inObsPtswHeight, heightfield, heightfield - heightOffset, "PYTHON")
TypeError: unsupported operand type(s) for -: 'unicode' and 'unicode'

EDIT: Bera have given around 70 % of the answer:
the first argument (which is a field) cannot be changed to float because it is a field name. the second argument however can.

heightfield = arcpy.GetParameterAsText(0)
#line above: field in table inObsPtswHeight - this is left unchanged
heightOffset = float((re.findall("\d+\.\d+",arcpy.GetParameterAsText(1)))[0])
# linear unit, constant value - a linear unit returns a string (e.g. "102.9 meters")
# which has to be stripped into only the number/float. this line returns a float.

arcpy.CalculateField_management(inObsPtswHeight, heightfield, "!{0}! - {1}".format(heightfield, heightOffset), "PYTHON")
# as Bera pointed out, the smartest is putting the expression into a string.
#to call on the field in a python field calculator expression, you have to use !fieldname!, while the float can be used directly.


texts: DXF export in QGIS giving hatches instead of polylines

texts:
I am new to QGIS, and using QGIS 2.0.1.

I'm having a problem when saving as .dxf. I export a building layer (polygons) as a .dxf file However in AutoCAD it turns up as a hatch instead of as a closed polyline which is what I'm looking for. I can of course recreate the boundary in AutoCAD, but I would rather just save the data correctly in QGIS.

Can anyone help?


texts: Can one batch-automate the creation of QGIS heatmaps?

texts: I have 4k point shapefiles, from which I would like to generate heatmap TIFFs using the QGIS heatmap plugin (or a functional equivalent featuring triweigth, radius and cell size selection). Is there any way to batch / script / automate this? I.e. "add vector shapefile, run heatmap plugin with specific options, remove vector shapefile and raster heatmap, add next vector shapefile, repeat..."?


texts: Is it possible to get column names from an ArcSDESQLExecute SELECT execution

texts:
I am running ArcGIS 10.0 desktop and querying an Oracle 11g database and writing the results of a SELECT statement out to a CSV file. Everything works great; however, I would love to find a way to include the column names in the CSV without hard coding them.

My code looks like this:

obj_sde_conn_ora = arcpyArcSDESQLExecute(server, port, "#", user, password)

sql_ora_select = """SELECT col1, col2
FROM tab"""

sde_return = obj_sde_conn_ora.execute(sql_ora_select)
csv_file = open(out_csv, "wb")
import csv
csv_writer = csv.writer(csv_file, delimiter=";", quotechar='"', quoting=csv.QUOTE_MINIMAL)
for row in sde_return:
    csv_writer.writerow([unicode(s).encode("utf8", "replace") for s in row])
del sde_return
csv_file.close()


texts: Is there PostGIS equivalent to arcpy.sa.Slice() function?

texts:
I have an existing system built on an 8.3 PostgreSQL database which uses a ArcGIS - Spatial Analyst for a couple of raster operations.  The system is moving to PostgreSQL 9.2 and I'm evaluating if the Spatial Analyst piece is actually needed.  But I'm not at all familiar with PostGIS's raster capabilities, and couldn't find an equivalent function during my initial search.

The actual function call is:
arcpy.sa.Slice(xras, zonecount, "NATURAL_BREAKS")

From the Esri Help Docs:

  Slices or reclassifies the range of values of the input cells into
  zones of equal interval, equal area, or by natural breaks.

  Slice works best on data that is normally distributed. When using    input raster data that is skewed the output result may not contain
  all of the classes that you had expected or specified.
  If an environment Mask has been set, those cells that have been masked will    receive NoData on the output slice raster.
  When using the EQUAL_AREA    method, sometimes not all of the output classes will have an equal,    or even similar, number of cells
  (i.e., area). This may be an    inherent result based on the nature of
  the input values and the    specified number of classes (bins). If the
  results are deemed    undesirable, you can try using a few number of
  classes or applying a    statistics transformation (e.g., logarithm,
  square root, and so on)    to the input dataset.


texts: Extracting centerline of a Complex-Polygon in PostGIS/Python

texts: Any idea how to extract the centerline of the following polygon (Fig1) using PostGIS functions or Python modules? I was trying Voronoi Function but thats not helpful for this kind of shape. I want something like Fig2 blackline.


texts:
I am building a modeler in QGIS 3.2.1, I want to use Refactor fields tool to change the name of a fixed field (fid to org_fid)，This tool accepts a Fields mapping. I also noticed that there is a Fields Mapper type parameter in the QGIS Parameters. If I add a Fields Mapper type parameter to the Modeler, I can specify this parameter in the Fields mapping of the tool. But I don't know how to specify fid to org_fid.
Can anyone provide an example of using Fields Mapper in the modeler?


texts: Running spatial analysis using ArcPy in PyCharm failed with error code 1073740791 (0xC0000409)

texts:
I'm running spatial analysis using ArcPy in PyCharm (Windows 10 64bit environment) failed with error code 1073740791 (0xC0000409). Specifically, i'm trying to run Expand within the spatial analysis module. My code is as followed:

import arcpy
import glob
import os
from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')

inpath=r'D:\Work\NBR\\'
outpath=r'D:\Work\NBR\expanded\\'

infilelist = glob.glob(inpath+'*.dat')
infile=infilelist[0] #using a single file as an example

basename = os.path.splitext(os.path.basename(infile))[0]
outfile=outpath+basename+'_expanded.tif'
if os.path.isfile(outfile) == False:
   zone_values = [-2500]

   outExpand = Expand(infile, 5, zone_values)
   outExpand.save(outfile)
   print outfile + ' is done!'

arcpy.CheckInExtension('Spatial')

I tested this script in the debug mode and it seems that every time when

outExpand.save(outfile)

is run, the program aborts and gives me a message:

  Process finished with exit code -1073740791 (0xC0000409)

No other error messages are given.

Can someone help me find out why this happens?

For your information, this same code can be run on my other machine without issues. I'm trying to get it to run on my newly-purchased Windows laptop. My GPU is an NVIDIA Quadro P1000. Its driver has been upgraded to the latest version 23.21.13.9103. The ArcMap version i'm using is 10.6.


texts:
I have a lot of .gz files that contain imagery and I'm going to unzip them.
I was wondering if that could be done via GDAL.


texts: Distance analysis of 3-D surface in ArcGIS using VB.NET

texts:
Using VB.NET I created a 3D polygon in ArcGIS to model a fault plane and now I want to find two different distances from a point on earth's surface. How to go about this? Can some one help with the code?

I have used following code to create polygon with 4 points. For example Lat1, long1, z1 are latitude, longitude and depth of the first vertex. Now I need to find the following distances:

Nearest distance from a point (Defined by its own lat and long) on earth's surface to the 3-D surface
Nearest distance from the same point to the polygon formed by projecting this 3-D polygon on earth's surface

Dim pPoint As IPoint = New PointClass()
Dim pPointCollection As IPointCollection = New RingClass()
Dim pGeomCollection As IGeometryCollection = New PolygonClass()
pPoint.PutCoords(lat1, long1)
pPoint.Z = Z1
pPointCollection.AddPoint(pPoint)
pPoint.PutCoords(lat2, long2)
pPoint.Z = Z2
pPointCollection.AddPoint(pPoint)
pPoint.PutCoords(lat3, long3)
pPoint.Z = Z3
pPointCollection.AddPoint(pPoint)
pPoint.PutCoords(lat4, long4)
pPoint.Z = Z4
pPointCollection.AddPoint(pPoint)
pGeomCollection.AddGeometry(pPointCollection)
Dim pOutGeometry As IGeometry = CType(pGeomCollection, IGeometry)
Dim pZaware As IZAware = CType(pOutGeometry, IZAware)
pZaware.ZAware = True

For the first distance, we can use "Near Tool" in ArcGIS but I am not knowing how to initiate that and find the distance.For the other distance we can form a 2D polygon by ignoring the Z values of vertices and use the same Near Tool.
Please help me with the code for this idea.


texts: How can I merge data back to SpatialPolygon (or retain initially) after nowrapRecenter command?

texts:
I am using newobj <- nowrapRecenter(obj), from the package maptools to recenter a SpatialPolygonsDataFrame such that longitude runs 0 - 360 (pacific view) rather than -180 to 180 (Atlantic view). Unfortunately, the resulting newobj is a SpatialPolygon only --- the data is lost. I expect I can merge the data back on, but I'm not sure how... generally I would merge by a data column, and this new polygon has no data.

I am using this GEOLEV1 shapefile:

https://international.ipums.org/international/gis_harmonized_1st.shtml


texts:
I had been working with BioClim data and that was working in Maxent. I could clip files and realign them etc.

Now I have some data from CliMond which is in two formats:

Either as text files

CM10_1975H_Bio01_V1.2.txt
CM10_1975H_Bio02_V1.2.txt
CM10_1975H_Bio03_V1.2.txt
etc

Or sets of data (lots of folders Bio1 to Bio40) as ESRI format files

/bio40$ ls

dblbnd.adf
hdr.adf
metadata.xml
prj.adf  sta.adf
w001001.adf
w001001x.adf

Is there a way that I can batch clip these to a format that is usable in Maxent (python or Qgis)?


texts:
I have an OpenLayers map and I try to implement a "lightbox" image gallery , in the same page.

By "lightbox" image gallery" I mean that the images/arrows/exit button cover the whole screen.

The problem is that OpenLayers map covers the lightbox gallery (see image).

I have set z-index to the map and to the lightbox divs, like so

<div id='map_element' style='width: 900px; height: 400px; z-index:1;'> </div>

//and then for lightbox
#imagegal2 img {
 height:85%;
 z-index:90;
}

And this does not work

How do I fix this?


texts:
I'm new to OpenLayers and looking for some help. I would like to filter my shapefile with CQL and display the result. (that's actually working :-)

The Problem is that I don't know how to pan the map to the result.
So here is my code. I would be overgrateful for your help!

suche = new OpenLayers.Layer.WMS(
    "nkGIS:HauptAdressen - Tiled", "http://localhost:8080/geoserver/nkGIS/wms",
    {
        LAYERS: 'nkGIS:HauptAdressen',
        STYLES: '',
        transparent: "true",
        format: 'image/png',
    },
    {
        buffer: 0,
        displayOutsideMaxExtent: true,
        isBaseLayer: false,
        yx : {'EPSG:900913' : false}
    }
);

map.addLayers([suche]);

osm = new OpenLayers.Layer.OSM();
map.addLayers([poi]);

map.addLayers([osm]);

poi.setVisibility(false);
suche.setVisibility(false);

// support GetFeatureInfo
map.events.register('click', map, function (e) {
    document.getElementById('nodelist').innerHTML = "Loading... please wait...";
    var params = {
        REQUEST: "GetFeatureInfo",
        EXCEPTIONS: "application/vnd.ogc.se_xml",
        BBOX: map.getExtent().toBBOX(),
        SERVICE: "WMS",
        INFO_FORMAT: 'text/html',
        QUERY_LAYERS: map.layers[0].params.LAYERS,
        FEATURE_COUNT: 50,
        Layers: 'nkGIS:HauptAdressen',
        WIDTH: map.size.w,
        HEIGHT: map.size.h,
        format: format,
        styles: map.layers[0].params.STYLES,
        srs: map.layers[0].params.SRS};

        // handle the wms 1.3 vs wms 1.1 madness
        if(map.layers[0].params.VERSION == "1.3.0") {
            params.version = "1.3.0";
            params.j = parseInt(e.xy.x);
                    params.i = parseInt(e.xy.y);
                } else {
                    params.version = "1.1.1";
                    params.x = parseInt(e.xy.x);
                    params.y = parseInt(e.xy.y);
                }

                // merge filters
                if(map.layers[0].params.CQL_FILTER != null) {
                    params.cql_filter = map.layers[0].params.CQL_FILTER;
                }

                OpenLayers.loadURL("http://localhost:8080/geoserver/nkGIS/wms", params, this, setHTML, setHTML);
                OpenLayers.Event.stop(e);

            });
        }

    function updateFilter(){
        if(pureCoverage)
            return;

        var filterType = "cql";
        var str = "SName LIKE'" + document.getElementById('str').value +"'" + "AND " + "HNR LIKE'" + document.getElementById('nr').value +"'";

        // by default, reset all filters
        var filterParams = {
            str: null,
            cql_filter: null,
            featureId: null
        };
        if (OpenLayers.String.trim(str) != "") {
            if (filterType == "cql")
                filterParams["cql_filter"] = str;
        }
        // merge the new filter definitions
        mergeNewParams(filterParams);
        suche.setVisibility(true)
    }

    function resetFilter() {
        if(pureCoverage)
          return;

        document.getElementById('str').value = "";
        updateFilter();
    }

    function mergeNewParams(params){
        suche.mergeNewParams(params);
    }


texts: Google Earth Engine Python API clipped export bad behavior

texts:
I am trying to export a multi-polygon clipped image with the python API. The problem is that the code in python does not clipped the image correctly (I have just one polygon saved instead of two like the example given below). Probably I am doing something wrong but I do not know where my mistake is. I applied the same code on the Google earth Engine editor with a similar Javascript code and everything is correctly done.

Javascript code:

 var geometry = ee.Geometry.Polygon([[[37.24154695,  8.01665693],
 [37.26876315,  8.01665693],
 [37.26876315,  8.04378808],
 [37.24154695 , 8.04378808],
 [37.24154695,  8.01665693]],
 [[37.34154695,  8.11665693],
 [37.36876315,  8.11665693],
 [37.36876315,  8.14378808],
 [37.34154695 , 8.14378808],
 [37.34154695,  8.11665693]]]);

/**
 * Function to mask clouds based on the pixel_qa band of Landsat 8 SR data.
 * @param {ee.Image} image input Landsat 8 SR image
 * @return {ee.Image} cloudmasked Landsat 8 image
 */
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterBounds(geometry)
                  .filterDate('2016-01-01', '2016-12-31')
                  .map(maskL8sr);
var image = dataset.median()
Export.image.toDrive({
  image: image.clip(geometry),
  folder: "test",
  description: 'imageToCOGeoTiffExample2',
  scale: 30,
  region: geometry,
  fileFormat: 'geotiff',
  skipEmptyTiles: true
});

Python code:

t = [[[37.24154695, 8.01665693],
          [37.26876315, 8.01665693],
          [37.26876315, 8.04378808],
          [37.24154695, 8.04378808],
          [37.24154695, 8.01665693]],
         [[37.34154695, 8.11665693],
          [37.36876315, 8.11665693],
          [37.36876315, 8.14378808],
          [37.34154695, 8.14378808],
          [37.34154695, 8.11665693]]]

    region = ee.Geometry.Polygon(t)
    dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
        .filterBounds(region) \
        .filterDate('2016-01-01', '2016-12-31') \
        .map(_mask_l8_sr)
    dataset = dataset.reduce('median')
    task = ee.batch.Export.image.toDrive(image=dataset.clip(region),
                                         description="test",
                                         folder="testing",
                                         region=region['coordinates'],
                                         scale=30,
                                         fileFormat='GeoTIFF',
                                         skipEmptyTiles=True)
    task.start()


texts: Selecting the transect lines which touches the divide of a polygon from the river line

texts:
I have a polygon and transect lines which passes inside the polygon.

Here I want to select the lines which touch the drainage and divide. The transect lines split at the river intersection.

Since all the transect lines lies inside the polygon the intersection selection selects most of the unwanted features. Is there any tool that I can use to automatically perform to get similar results?


texts: Srver-side (WMS) rendering / caching of vector data

texts:
I have to develop a webapp that will have to display a large number of polygons on a webpage, with complex filtering options. The source data never changes (I will simply reload a new version of the DB every few months).
I will probably use OpenLayers on the client-side. I developed a prototype that drawed the feature as vectors in javascript, but I'd like to switch to server-side rendering for better performance. I could probably do my app connecting directly OpenLayers to a few different WMS layers in GeoServer (postgis backend) with CQL filters, but I have the feeling GeoServer is not very flexible when using a combination of complex filters and dynamic styling of my WMS and the solution could become quite complex.
From the developer point of view, a more confortable/flexible approach seems to use GeoDjango and its spatially-enabled data models to generate a dynamic (content and styling) KML file of my data. A (probably naive) elegant solution would be to integrate WMS rendering (and possibly caching) layer between the OpenLayers client and the GeoDjango app that generates KML output.
The workflow would be:

Openlayers sends user query to "intermediate layer"

This layers forward query to the GeoDjango app, get a KML file and render it as WMS and return it to OpenLayers.

(alternatively, OL could also talk directly to GeoDjango it sould use a middleware to postprocess its KML result before returning it).
As my data almost never changes, I could then benefit from a caching tool such as GeoWebCache between OL and the "rendering layer".
I know Google Maps API do server-side rendering of KML files when they are large, but it's a black box tied to their JS library, and I'd really like to do something similar with OL.
I've done a few searches, but was not able to find any tool for the "intermediate layer". WMS servers such as GeoServer/Mapserver can use a KML data source, but it seems limited to a static file somewhere in the filesystem, and not from a dynamic URL.
Does a tool like this exists or is my approach simply wrong ?


texts: Why do polylines digitised disapper when I save the layer. Would the shapefiles be corrupt?

texts:
I have a problem with QGIS 2.2. When I digitise a polyline layer and save the edits the lines disappear from the map. This does not happen for the polygons in the same GIS project. The points are still present in the attribute table.

Also related problem with the lines which have been categroised to show different symobology do not display with the different symbology

I have never had this problem before. Could the files be corrupt?


texts: How can we calculate the area of dense and sparse vegetation from NDVI image in QGIS?

texts: I'm trying to obtain values of area covered by dense and sparse vegetation through NDVI classification. I'm using QGIS (2.18), ERDAS 2014 and ArcMap.


texts: QGIS change zones in a raster to grayscale by polygon overlay

texts:
I've got a raster ortophoto and a polygon layer:

I need to alter the raster layer as follows:

Everything that is inside the polygon (overlaying) should stay as it is right now (in color).
The rest of the raster image (the zone that is not overlaying the polygon) should be changed to grayscale.

I can't imagine right now how to achieve this..


texts: Split area into equal size squares, count number of records for each segment

texts:
I'm working on a project where I have a series of points, something like this:

| id   | name        | geom  |
|------|-------------|-------|
| 1    | example a   | point |
| 2    | example b   | point |
| ..
| 641  | example bcc | point |

Given a bounding area (always rectangular), such as:

{
    "type": "Polygon",
    "coordinates": [
      [
        [
          -0.23895263671874997,
          51.390637019843744
        ],
        [
          0.06866455078125,
          51.390637019843744
        ],
        [
          0.06866455078125,
          51.59413527860747
        ],
        [
          -0.23895263671874997,
          51.59413527860747
        ],
        [
          -0.23895263671874997,
          51.390637019843744
        ]
      ]
    ]
  }

I wish to divide the area into 100 equal size segments, and calculate the total count of records for that area (the intention being when a user zooms in, the bounding box can be narrowed and more specific results provided).

For example, returning something like this:

| geom    | count |
|---------|-------|
| polygon | 26    |
| polygon | 12    |
| polygon | 8     |

I haven't been able to find any function for splitting a polygon other than by a Line, however this would require creating 18 lines (9 vertical, 9 horizontal), then looping over them. Is there a better way?


texts:
I have an irregular polygon.

I want to make parallel lines with 14 meters of space between the inside of the polygon using QGIS, does not matter the direction, but it is necessary to have the least lines as possible.


texts: MapInfo - programmatically find if layer is hidden

texts:
I am working on an integrated mapping project. I would like a button that toggles a layer. I could hold the layers status in a variable and then show/hide it as appropriate but this seems a little hacky.

Is there a way to find out whether a layer's "display" value is set to on or off?


texts: How to solve the problem when buildings are extruded without roof?

texts:
I tried to extrude buildings in ArcScene with Extrusion function in Layer Properties. But most of the buildings seem like only their walls are extruded. Instead of solid blocks, they appear like no roof and ground.

Does anyone know how to solve this problem?


texts: Update raster data without restarting the service

texts:
I am trying to update raster datasets (3-4 times a day) that are being used by an ImageServer service.

Is there a way to do this without stopping/restarting the service (to get rid of possible locks) ?

PS: I am using ArcGIS Server 10.3.1 and the files are .tiff.


texts: How to calculate mean of one band of an image in rgee (R package for interacting with GEE) without exceeding memory limits?

texts:
I am calculating the means of an image that I've previously done calculations on but I am exceeding EE memory limits: "Error in py_call_impl(callable, dots$args, dots$keywords) :
EEException: User memory limit exceeded."
Is there a more efficient way to calculate the mean?
means1 = image$reduceRegion(
  reducer=ee$Reducer$mean(),
  geometry=shp$geometry(),
  scale= 10,
  maxPixels= 1e9)
print(means1$getInfo())


texts: Using pgr_withPointsCost from pgrouting 2.2.2

texts:
I'm trying to calculate the cost of travelling in a network from one point to many points (adresses).
To increase the precision of the analysis i want to use the function pgr_withPointsCost as it (as far as I understand) can calculate to and from points that are not the actual vertexes of the network.

I have loaded osm data into postgres using osm2pgrouting and I have created a table containing the id's of my points of interest (poi), the id of the closest edge (edge_id) and the projected location of my poi on the edge as a fraction.

When running the pgr_withPointsCost query I get an output of 0 rows. Below is my create table statements and the current query - can anyone help me find the problem with my query?

ways table from OSM2pgrouting :

CREATE TABLE osm_network.herlev_ways
(
  gid bigint,
  class_id integer,
  length double precision,
  length_m double precision,
  name text,
  source bigint,
  target bigint,
  x1 double precision,
  y1 double precision,
  x2 double precision,
  y2 double precision,
  cost double precision,
  reverse_cost double precision,
  cost_s double precision,
  reverse_cost_s double precision,
  rule text,
  one_way integer,
  maxspeed_forward integer,
  maxspeed_backward integer,
  osm_id bigint,
  source_osm bigint,
  target_osm bigint,
  priority double precision,
  the_geom geometry(LineString,4326)
)

Points of interest table:

 CREATE TABLE base_data.geo_borger_aktiv
(
  ogc_fid serial NOT NULL,
  wkb_geometry geometry(Point,25832),
  kvh_adr_key character varying,
  CONSTRAINT geo_borger_aktiv_pkey PRIMARY KEY (ogc_fid)
)

Nearest network edge to poi - including fraction:

CREATE TABLE base_data.nearest_ways_borger_aktiv
(
  id bigint,
  pid character varying,
  edge_id bigint,
  fraction double precision
)

Query:

SELECT *
FROM pgr_withPointsCost('SELECT gid AS id, source,
target, length_m AS cost, length_m AS reverse_cost
FROM osm_network.herlev_ways
WHERE class_id > 105
ORDER BY gid',
'SELECT id*-1 AS pid, edge_id, fraction FROM
base_data.nearest_ways_borger_aktiv',
117558
,
(SELECT array_agg(nearest_ways_borger_aktiv.id*-1)
 FROM base_data.nearest_ways_borger_aktiv
 WHERE nearest_ways_borger_aktiv.id IS NOT NULL))

UPDATED withPointsCost query:

SELECT *
FROM pgr_withPointsCost('SELECT gid AS id, source, target, length_m AS cost, length_m AS reverse_cost
FROM osm_network.herlev_ways
ORDER BY gid',
'SELECT id AS pid, edge_id, fraction FROM base_data.nearest_ways_borger_aktiv',
-150980,
(SELECT array_agg(nearest_ways_borger_aktiv.id*-1)
 FROM base_data.nearest_ways_borger_aktiv
 WHERE fraction BETWEEN 0.0001 AND 0.9999
 ))

Creating edge_id and fraction:

CREATE TABLE base_data.nearest_ways_for_busstops AS

SELECT DISTINCT ON(ogc_fid)
ogc_fid AS id,
gid AS edge_id,
ST_LineLocatePoint(ver.the_geom,ST_Transform(akt.wkb_geometry,4326)) AS fraction
FROM osm_network.herlev_ways ver, base_data.nearest_ways_borger_aktiv akt
WHERE ST_DWithin(ver.the_geom, st_transform(akt.wkb_geometry,4326), 0.001)
ORDER BY ogc_fid, ST_Distance(ver.the_geom, ST_Transform(akt.wkb_geometry,4326)) ASC

Points of interest (nearest_ways_borger_aktiv) visualised, including edge_id of nearest ways and calculated fraction, along with ways with corresponding gid:


texts: Using Add Field, Calculate Field within ModelBuilder?

texts:
I have two questions regarding ModelBuilder.

I have a point layer with Crop Yield values. I need to Add Field called "Count", then do a Calculate Field, Count = 1. The problem that I'm running into is that in Model builder, when I try to insert a Calculate Field tool, "Count" field does not exist yet so I cannot reference that. How do I get the Count field to appear in MB before running the Calculate Field operation?
I need to do a Spatial Join in ModelBuilder as well. I have done a Multiple Buffer operation with 7 concentric circles, and I need to join each field (each "doughnut) in that layer with the points that fall within it. I then need to perform several operations including Sum, Min, Max, Std Dev. When I manually do a "right-click" Spatial Join (right click the layer itself and select Join), this is no problem. Here's the dialog I would like to see in Model builder:

However, the Spatial Join tool that is brought into ModelBuilder is significantly different. In it, you have the list of fields from the Join Feature, and you have to select an individual operation for each one.

Can someone tell me how to do this Join in ModelBuilder?


texts: Linear referencing with different "begin/end" offset

texts:
Here is my problem: I have an axis passing in the middle of the street that would be used for linear referencing.

Now I need to draw each border of the road.

On the picture, I have no problem with the right part since the width of A and B is exactly the same, but is it possible to adapt the offset for the beginning and the end of the line to represent left side accordingly to the picture ?


texts:
I have used the line intersection tool in QGIS 3 and I have used vector id for example: vector 1 intersects with 3 and the point (id 1) but also repeats the same imformation in reverse saying vector 3 intersects with vector 1 and then gives another point (id 2).

Is there a way of deleting or sorting the attribute table so the information is not repeated?


texts:
using QGIS 2.18.1 in Windows 10, what I have is 2 raster files with different cell sizes.
raster1 has small cell sizes, and represents land uses (each cell value represent a land Use) and raster2 has bigger cell sizes and represents precipitation values. You can see it from the image:

I need to assign to each cell in raster1, the value of the overlapping cell value of raster2. And in the cases in which a cell in raster1 overlaps with 2 cells in raster2 (A cell between 2 cells of raster2), get the value of the raster2 cell overlapping bigger surface.

I could do this by polygonizing both raters and applying python code to get most predominant values within polygons (like in this question) but I wonder if there's an easier way to do this by using the raster calculator or other raster tools.


texts:
QGIS 2.18. (Las Palmas)

I would like to create a series of maps, where I rather than zooming to the coverage Layer, i would only highlight the filter i used (for example with a different color) but still want to be able to see the whole map. Is there a possibility to do that?


texts: ArcGIS for Server 10.1, Registering a database data store fails with Underlying DBMS error

texts:
When I registered a database data store to ArcGIS Server in ArcCatalog, the following error occurs:

  The connection property set was missing a required property or the
  property value was unrecognized. Underlying DBMS error [Unknown
  errorNo extended error.]

I set up the ArcGIS Server and Oracle on two computers as follow

192.168.1.51

ArcGIS for Server 10.1 64bit
ArcCatalog 32bit
Oracle Instant Client 11.1.0.7 32bit
Oracle Instant Client 11.1.0.7 64bit

192.168.1.52

Oracle Database 11g Enterprise Edition 11.1.0.7.0 64bit

I have successfully made the Oracle 11g "Enabled Enterprise Geodatabase" through the ArcCatalog, and it seems nothing wrong.

I added the publisher database connection with the following variables

Database Platfrom: Oracle
Instance: 192.168.1.52:1521/orcl
Authentication Type: Database authentication
User name: sde
password: sde
[v] Save user name and password

And I have added the Oracle Instant Client directory to the PATH and environment variable, I am sure that the 64-bit DLL is found before the 32-bit one. But I still get the error message.

I am wondering if there is anything else I have to set up?


texts: How to upload a large amount of datasets into GeoNode?

texts:
I'm wondering to upload a large amount of datasets into GeoNode (from shapefiles or PostGIS)... but if I add those datasets directly to GeoServer, GeoNode won't know that a specific layer was added and for consequent I wont see in GeoNode front end my uploaded layers..

How is the best way to do this?
Design a script that uses GeoNode classes to upload shapefiles? And what about if I want to link a database directly to GeoServer for automatic layers creation? Is there any way to GeoNode look for updates that were made in GeoServer and show them on GeoNode front end?

It must be develop a new Django application to do this?

What is the best approach?


texts: I have a table that contains buildings in a line geometry where a line starts and ends at the same point. Is there away to convert lines to polygons in Oracle Spatial?


texts: ECW file solution doesn't work in QGIS in the latest ubuntu

texts: In Ubuntu 18.4 I enrolled in a QGIS course in that course .dem file shows but .ecw can not show so what can be the solution can anybody tell me. QGIS version is 3.10.


texts:
I have a geotiff file and I use geotools library for parsing it.

This is my code:

File f = new File("myFile.tif");

ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
policy.setValue(OverviewPolicy.IGNORE);

ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
useJaiRead.setValue(true);

GridCoverage2D image = new GeoTiffReader(f).read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});

Rectangle2D bounds2D = image.getEnvelope2D().getBounds2D();
GridGeometry2D geometry = image.getGridGeometry();

System.out.println(geometry);

System.out.println(bounds2D);

The output of geometry is:

GridGeometry2D[
    GeneralGridEnvelope[0..1199, 0..1399],
    PARAM_MT[
        "Affine",
        PARAMETER["num_row", 3],
        PARAMETER["num_col", 3],
        PARAMETER["elt_0_0", 996.4799194335938],
        PARAMETER["elt_0_2", -598386.1975402832],
        PARAMETER["elt_1_1", -999.6777954101562],
        PARAMETER["elt_1_2", 650290.4111022949]
    ]
]

The output of bounds2D is:

java.awt.geom.Rectangle2D$Double[x=-598884.4375, y=-748758.6635742188, w=1195775.9033203125, h=1399548.9135742188]

What does this information mean?


texts: Creating curved lines between two points in QGIS

texts:
I am trying to create curved lines between two point in QGIS. There are a few posts that show how to create these lines (hyperlinks below). However, I am having trouble executing the code in QGIS.

This article says one could use the following code:

SELECT st_makeline(st_segmentize((SELECT geog FROM paths WHERE gid = 1), 50000)::GEOMETRY);

This article says one could use this code:

UPDATE experimental.airroutes
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
       ST_Transform(a.the_geom, 953027),
       ST_Transform(b.the_geom, 953027)
     ), 100000 ), 4326 )
 FROM experimental.airports a, experimental.airports b
 WHERE a.id = airroutes.source_id
   AND b.id = airroutes.dest_id
);

It's unclear to where where to execute this code. I have tried several options that I outlined below which so far have not been successful. Here is what I did so far:

I think I may need to utilize the "Add PostGIS Table(s)" (the elephant icon) function. If so it's unclear to me what the Name, Service, Host, Database, Username and Passwords should be. This post asked about this and the respondent suggested accessing pgAdmin, which I do not have and can't figure out how to get.

This post suggests executing the PostGis SQL code from the PostGIS Manager Plugin. However, this is not an available plugin from the Plugins library.

I have also tried to use DB Manager but when I click on the SQL Window icon I receive a message noting "No database selected or you are not connected to it."

Finally, I installed the Post GIS SQL Editor plugin but when I open the plugin I receive the following error "An error has occured while executing Python code:

TypeError: QComboBox.findText(QString, Qt.MatchFlags flags=Qt.MatchExactly|Qt.MatchCaseSensitive): argument 1 has unexpected type 'QPyNullVariant'


texts: Merge two point datasets and selecting the greatest / smallest Z value

texts:
Is there a way to merge two point value datasets which overlap, with an argument to always take either the lowest or highest Z value at points which overlap? I want to create a hybrid dataset of two overlapping layers, with the argument to select the greatest value for Z where they overlap.

I understand that Cell Statistics can be used for this, but that requires raster inputs and outputs, I need the data to remain in an XYZ point dataset format so that it can then be exported to another software package that does not handle rasters.

I can merge the two point datasets into a shapefile using the Merge function, however this has no option to add an argument for overlapping sections.


texts: Converting list of polygons to multipolygon using shapely?


texts: Is it possible to get arcpy on 9.3* or Qgis? Or maybe some of the functions?

texts: If this isn't possible what would be the closest equivalent open source? Could parts of the module be ported to work on 9.3 or QGIS? I'm not sure if this would infringe on the license porting for 9.3.


texts: How to avoid "Unable to compute a transformation between pixel/line"?

texts:
I am trying to re-project from Equi-rectangular to the NASA global EASE grid (Equal Area Scalable Earth) grid 1383 586 with a ~25 km cylindrical projection (ESPG=3410). The file is provided as HDF5.

But I am facing an error:

  gdal_translate -of "ENVI" HDF5:"C:\Users\GW1AM2_20120930_01D_EQMD_L3SG
  SMCHA1110100.h5"://Time_Information C:\Users\amsr2_time.img

  gdalwarp  -of "ENVI" -overwrite -ts 1383 586 -t_srs "+proj=cea +lon_0=0 +lat_ts=30 +x_0
  =0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs" -ot Float32 C:\Users\amsr2_time.img C:\Users\amsr2_time_ease.img

ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates for C:\Users\amsr2_time.img.
There is no affine transformation and no GCPs.

The file can be downloaded as both HDF5 and NETCDF` formats here:

http://gcom-w1.jaxa.jp/OrderFiles/31DaKuDQzFDeAgM0/31DaKuDQzFDeAgM0_0.zip


texts: Creating raster/DEM that represents pre-landslide morphology using ArcGIS Desktop

texts:
I've posted something similar here before, but my previous post did not yield great results.

I am a Master student in Geology studying landslides. For a part of my project I need to calculate the volume of a several prehistoric landslide deposits. Most of the papers I have found so far discuss landslides deposits which have occurred relatively recently and have the benefit of some kind of pre-slide topography, satellite image, interpolated GPR, or some other kind of DEM, making the cut-and-fill method (ArcMap) for calculating volume fairly straight forward. Since the slides I'm looking at are all at least 1000 years old, there is no satellite imagery or topo. map to pull from.

Fortunately I do have access to an excellent Lidar dataset in my study area for the post-event DEM, but my issues arise when trying to construct a pre-event raster/DEM. Specially I'm looking for a step-by-step guide on how to do this (my ArcGIS skills are subpar at best) so I can insert my custom pre-event DEM into the cut-and-fill tool. I understand that I essentially just need to create a DEM that represents what the topography probably looked like before the landslide; I just don't know how to do this for ArcGIS.

I will have to calculate both the volume of the debris field and the volume of the source zone (depletion zone).

(1) Volume of debris field

Since these landslides are all located in a river valley, there are river terrace deposits which are relatively flat. The terraces can be used to get an idea of the pre-event topography. They appear well on LiDAR but my problem is accounting for small elevation changes from terraces further upstream and those downstream. I am not sure how to create a slightly tilted DEM that would represent a pre-event river terrace.

(2) Volume of source zone

Second is calculating the void space or the amount of material that was released during the landslide event. The hollows are fairly evident when looking at LiDAR, and theoretically I would just "drape" a DEM over the void space like holding a sheet taut over the hollow.

The best paper I have found so far (https://www.nat-hazards-earth-syst-sci.net/18/709/2018/) mentions "the depleted area was manually filled by inserting several linearly aligned points within the slide area" (section 4.2 in results). I just don't know enough about ArcGIS to figure this part out.


texts: Addressing "Invalid Manifest" error when sharing geoprocessing tool script to ArcGIS Enterprise as a web tool

texts:
I wrote a Python script which takes an input Cropland Data Layer raster dataset and a user's defined Region of Interest, and returns a pop-up message telling the user the number of acres of each crop in that ROI. The script runs fine in ArcGIS Pro 2.5 as a script tool.

Ultimately, I want to now share this script tool as a web tool on enterprise. However, when I try to do this, the process fails out, and when I check the logs on the server, the message which seems to indicate the failure is

  Invalid manifest. Unable to find dataset in extracted geodatabase.

I'm new to using Enterprise and Server and don't know how to address this error (I'm also not a Python wizard! Please forgive any poor coding). I would love any feedback pointing me in the right direction as to how to address this error so I can get the tool shared as a web tool. The code I have written follows.

I have a feeling the error lies in how I'm passing intermediate outputs, but I would not be surprised if I'm wrong.

Edit: I forgot to mention that I can successfully share the script tool as a geoprocessing package, if that matters

import arcpy
from arcpy.sa import ZonalHistogram
import tkinter
from tkinter import messagebox

#Set environments
arcpy.env.workspace = r"C:\Users\XXXXXXXXX\Documents\ArcGIS\Projects\CropAcreageTool\CropAcreageTool.gdb"
arcpy.env.overwriteOutput = True
root = tkinter.Tk()
root.wm_withdraw()

#Begin by having the user input their two parameters: Raster layer and ROI

#Parameter 1: Raster layer (input)
baseRaster = arcpy.GetParameterAsText(0)

#Parameter 2: Region of Interest (input)
ROI = arcpy.GetParameterAsText(1)

#Intermediate: Crop acreage table
cropTable = r"cropCount"

#Class key for join
classTable = r"classKey.csv"

#Execute ZonalHistogram to give the counts for each unique class within the ROI in the output table
ZonalHistogram(ROI, "OBJECTID", baseRaster, cropTable)

#Join the output table with the class key table to associate counts with crop type
arcpy.JoinField_management(cropTable, "OID", classTable, "OID")

#Loop through table with a search cursor and perform calculation to generate acreage: (count * 900m2)/4047m2 per acre
#Write result to list
#Display list as a popup message
cursor = arcpy.da.SearchCursor(cropTable, ["OBJEC_1", "Class_Name"], """"OBJEC_1" > 0""")

Results = []

for i in cursor:
    Results.append("{} acres of {}".format(round((i[0] * 900) / 4047), i[1]))

messagebox.showinfo("Results", "\n".join(Results))


texts: Creating an Interactive Chart/Graph in ArcGIS Online

texts:
My question: Can this interactive chart element be done in ArcGIS Online with some clever CSS additions? Or would I need to host the map with ArcGIS for Server and use my own JavaScript?

In ESRI's The ArcGIS Book, they show a dynamic map of tornado occurance as an example of how ArcGIS Online (AGOL) can be used. This map is available here. When looking at the "behind the scenes" of this product (here), they list ArcGIS Online as the required software--however AGOL doesn't allow for custom JavaScript. You would need ArcGIS for Server to host it.

In the side panel, they have an interactive chart. When a chart element is clicked, the map reacts (see image below). The chart summary data also updates.

My question: Can this interactive chart element be done in ArcGIS Online with some clever CSS additions? Or would I need to host the map with ArcGIS for Server and use my own JavaScript?


texts: My projection coordinate system does not show right figures on my print layout axes

texts: I am using QGIS 3.10.2 and trying to figure out why my coordinate system, although it shows right Lat/long figures on basic working screen(see 1st photo), when I try to work on Print Layout gives me wrong Lat/Long values. Small island on map(red circle) should be 16'5739.55 43'96496.18 under Greekgrid GGRS87 EPSG:2100 or 22.38181.886268 48.16018.276529 under EPSG:3857 WGS84 instead it shows 10'00 N and 16'E. All project layers are setup to EPSG:2100. All WMS layers are correctly aligned with my project.Project CRS is set to EPSG:2100. Appearence CRS is set to EPSG:2100. In PrintLayout Map/main propereties/grid is also set to EPSG:2100.


texts: Shortcut for zoom to data driven page extent?

texts:
I would like a shortcut to return my (data) view to the current data driven page extent.  Somewhat like having a bookmarked view, if the bookmark updated with the current data driven page.

Since I don't think such a thing exists, I thought I would create a customized keyboard shortcut.  But I haven't any idea what command I could set the shortcut for.

My work-around, and the behavior I would like to emulate, is to put the cursor in the text box on the Data Driven Pages toolbar, and then move it elsewhere - click somewhere else or press enter.  Anyone know how to put this into a keyboard shortcut?

Perhaps I should add that I'm using ArcMap 10.4.1.


texts:
In Postgres with PostGIS extension, I try to figure out to find the path containing a specific point:

table A contains POINT
table B contains PATH

My code returns 0 rows when it should

select * from A join B on ST_Contains(A.coordinates::geometry, B.coordinates::geometry);

I also tried to inverse A in B, without success. For the both tables ST_IsValid(coordinate::geometry) return true


texts:
I'm running QGIS 2.6 on a Mac, and trying to create buffers around a point shapefile. The process seems to work: (Vector>Geoprocessing Tools> Buffer(s), I select my input layer, leave the segments at 5, and set the distance, create a new output shapefile in the same folder as the input layer, and set the new layer to appear.

A layer IS created, but it's an ovoid ellipse that is massively too large - when I set it for '500 metres', it came out covering most of Eurasia and Africa! I tested the settings, and to at '0.01 metres', it produced an ellipse measuring approx 1200x1450 metres.

What is causing this, and how do I fix it?
I wonder if it's related to the co-ordinate system I'm using, as the project is WGS84-PseudoMercator and the input shapefile is WGS84 (I didn't intentionally make them different).

Any help gratefully received!
Adam


texts: How to extract latitude/longitude values for specific points from a shapefile?

texts:
I have obtained public crash data for the state of Minnesota.  I would like to be able to programmatically extract the latitude/longitude of each crash using this data in combination with a series of shapefiles that cover all roads in the state.  For each crash, I believe I have enough information to do this, but my knowledge of GIS is keeping me from getting things done.

I'll start with a simple example in the hopes that someone can point me in the right direction.  You'll have to forgive me if my GIS terminology is not what it should be.  The shapefile I'm using is here:  http://www.dot.state.mn.us/maps/gisbase/datafiles/county/ramsey.zip

Accident 090710194, for example, should be at mile 1.228 along TIS_CODE 0462000015.  When I use a tool like Quantum GIS to find all attributes with a TISCODE of 0462000015, there are 81 features to choose from.  If I select the feature with a BEGM (beginning mile) and ENDM (ending mile) that encapsulates 1.228, I find one feature that goes from BEGM of 1.21 to ENDM of 1.3.  This puts me "in the ballpark" for a precise location.  Then, I can use Quantum GIS' "Coordinate Capture" plugin, to guess that it's somewhere around -93.20228,45.05062.  I know these measurements are right, because a double check on Google Maps shows the same approximate location as Quantum GIS:  http://maps.google.com/maps?q=45.05062,+-93.20228

Really, there are 2 problems with the above approach:  1) It's very slow; I need to do this for over 500K crashes, so I need to do something programmatically, and 2) It's not very precise.

I've been told by folks in the Minnesota Department of Public Safety that if I had ESRI ArcMap (the tool they use), I could "create a point event", which sounds like a precise measurement along the feature I selected, and from there I would see the precise location, and (I'm hoping) could probably get latitude/longitude.  I don't have this software, but it sounds something like what is shown in the following video, where toward the beginning of this video they're creating a "route event".  Since I'm not allowed to post more than 2 links in this post, in order to see what I mean you'll need to Google "arcgis point event" and select the first link, which is to webhelp dot esri dot com.

Everything boils down to 2 questions for me right now:  1) Is programmatically extracting this information possible, even with ESRI ArcMap, and 2) If so, what (preferably open source) tools can aid me?  PostGIS?  Geotools?  Geoserver?  Any advice at all would be appreciated.


texts: Conditionally coloring specific country using Leaflet?

texts:
I have an offline Leaflet map, where the countries are drawn from a huge Javascript variable with coordinates.

index.html is as such(an example from codepen, I load the libraries locally):

  <title>Leaflet</title>

  <meta charset="utf-8" />
  <meta name="viewport" content="width=device-width, initial-scale=1.0">

  <link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />

  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.3/dist/leaflet.css" />

And Map.js is as such:

var map = L.map('mapid').setView([0, 0], 2);

// Huge geojson of every country, one country sample here.
var countries =  "type": "FeatureCollection",
"features": [{
    "type": "Feature",
    "id": "AFG",
    "properties": {
        "name": "Afghanistan"
    },
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [61.210817, 35.650072],
                [62.230651, 35.270664],
                [62.984662, 35.404041],
                [63.193538, 35.857166],
                [63.982896, 36.007957],
                [64.546479, 36.312073],
                [64.746105, 37.111818],
                [65.588948, 37.305217],
                [65.745631, 37.661164],
                [66.217385, 37.39379],
                [66.518607, 37.362784],
                [67.075782, 37.356144]
                // .... goes on
            ]
        ]
    }
}]};

L.geoJSON(countries, {
  style: function(feature) {
    return {
      fillColor: "#D3D3D3", // Default color of countries.
      fillOpacity: 1,
      stroke: true,
      color: "grey", // Lines in between countries.
      weight: 2
    };
  }
}).bindPopup(function(layer) {
  return layer.feature.properties.name;
}).addTo(map);

This produces a nice looking output, which can be seen here: http://codepen.io/anon/pen/zZNjBy

However, I wish to color potentially each country depending on an input. Let's say we open an object with the potential input:

let colorratio = {
    1:"green",
    2:"blue",
    3:"yellow",
    4:"red"
};

I then receive the following:

let colorinput = {
    AFG: "1",
    DNK: "2",
    SWE: "3",
    USA: "1"
};

Where we in layman's terms want to color Afghanistan green, color Denmark blue, Sweden yellow, and USA blue. The remaining countries should just stay at the default color value(greyish).

These abbreviations match the "ID" key in the countries variable.

My question therefore goes on how to approach this issue; Obviously, I need to match the keys in the huge GeoJSON file. But it's almost 13000 lines long.

Apart from that, in the L.geoJSON where we call the countries variable, we return a default color of countries. Where do I potentially access a key to color each country? Should this variable be moved into the countries variable for each country, as a key that is then overwritten?

I've tried looking at the Leaflet documentation and examples; but most assume that you retrieve tiles online, which is not a possibility here. It must work offline without contacting a CDN or similar - which poses some challenges, needless to say.


texts: Placing a label above a feature in 3D/Scene using ArcGIS Pro

texts: Is it possible (using ArcGIS Pro, 2.5.1) to manually or by derived extrusion (i.e. from a DTM) set the height of a label to sit/float above a feature? I see you can use an extrusion but that requires the label vector having a z value?


texts: What simplifying polygons functions f y = f(x) have ( y union x equal y ) property?

texts: I've looked at GRASS v.generalize manual but they don't say. As far I can tell there's hermit but it's a Smoothing function.


texts: Help with trying to track ice and polar bears on a map. No idea where to begin

texts:
I was asked create a little google map showing the location of Polar bears. I've been given the location data in lat-long format.

For each bear I've been given a file containing the following headings:

Date    Time    Year    Julianday   Hour    Minute  Activity    Temperature Latitude    Longitude

The intention was to show the migration of these bears over time - so I'd take the date, take the time, take the lat / long and animate a google map (straight ahead javascript  + googlemaps api stuff - got it working... easy).

However, I've now been asked to show the movement of arctic ice as well as bears.... Google maps doesn't have a North Pole... and while I've managed to overlay some ice, it just looks like a big white blob which didn't render out properly due to the lack of North Pole.

I've been asked to use data found here http://nsidc.org/data/masie/.

I need to find a way to take the nsidc data, and take the polar bear data and combine them to form a series of png files.

Once I have the png files I can animate them.... vistor hits 'play' and all the files play back to show the bears and ice moving over time. The play stuff is easy... I just don't know how I am to create these files...

To compound things, I've also been asked to include shipping data (also supplied in lat lon format)

Normally I'd do lots and lots of reading to try and figure stuff out- however, I've been given 3 days to do this... I've tried a few things...downloaded "QGis" and pulling in the data - but nothing appears.  I don't have a clue.

How do I / where do I begin?

Sorry - this is the worst, most vague question ever....


texts: ArcObjects ExcelWorkspaceFactory Won't Release Excel File

texts:
I'm working with ArcObjects to set up a batch import from xls and csv files into a geodatabase. That part works fine, but after the import, I want to delete the file. This, however, doesn't work, because the file is still open. None of the interfaces being used descend from IDisposable or expose a Close method. I've also tried using a ComReleaser to force close all of the COM objects that I open, but that doesn't work either. The only things I've seen pertaining to closing a connection say that I need to remove all references to these types and force garbage collection, which is against standards (for good reason).

Here's the method I'm using for the excel import:

private ServiceResult<IEnumerable<int>> InsertControlsFromXls(FileInfo xlsFile) {
  ServiceResult<IEnumerable<int>> result = new ServiceResult<IEnumerable<int>>();

  using(ComReleaser comReleaser = new ComReleaser()) {

    //read in the spreadsheet as a workspace
    Type factoryType = Type.GetTypeFromProgID("esriDataSourcesOleDB.ExcelWorkspaceFactory");
    IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);

    IWorkspace workspace = workspaceFactory.OpenFromFile(xlsFile.FullName, 0);
    IFeatureWorkspace remoteWorkspace = (IFeatureWorkspace)OpenWorkspace();
    IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
    IWorkspaceEdit2 session = (IWorkspaceEdit2)remoteWorkspace;
    IMultiuserWorkspaceEdit multiuserWorkspace = (IMultiuserWorkspaceEdit)remoteWorkspace;

    comReleaser.ManageLifetime(workspaceFactory);
    comReleaser.ManageLifetime(workspace);
    comReleaser.ManageLifetime(remoteWorkspace);
    comReleaser.ManageLifetime(featureWorkspace);
    comReleaser.ManageLifetime(session);
    comReleaser.ManageLifetime(multiuserWorkspace);

    try {
      //session.StartEditing(false);
      multiuserWorkspace.StartMultiuserEditing(esriMultiuserEditSessionMode.esriMESMVersioned);
      session.StartEditOperation();

      Match tableMatch = _tableRecReg.Match(xlsFile.Name);
      string tableName = tableMatch.Success ? tableMatch.Groups[1].Value + "$" : "Sheet1$";

      ITable table = featureWorkspace.OpenTable(tableName);

      result.Value = CopyAllFeatures(table, remoteWorkspace.OpenFeatureClass(AppSettings.DefaultSchema + "." + ControlMap.TableName));

      session.StopEditOperation();
      session.StopEditing(true);

    } catch {
      session.AbortEditOperation();
      session.StopEditing(false);
      throw;
    }
  }

  return result;
}

I asked this question on Stack Exchange too, but no responses so far:
https://stackoverflow.com/questions/30783703/arcobjects-excelworkspacefactory-wont-release-excel-file

Edit: To clarify, this is using ArcObjects 10.3, for which the API is slightly different from 9.3.


texts: In Google Earth Engine, I let the user define a polygon to give the geometry that the routine is applied to. As part of that, I would like to determine a "center" of the polygon. How does one extract a center of an irregular polygon?


texts: QGIS: error in conversion from vector to raster

texts:
I've having troubles with the raster conversion from a vector polygonal file.

This is the error in the log window:

INFO    GDAL execution console output
             ERROR 1: Wrong value for -outsize parameter.

Does someone have any advice?


texts: Why is spatialite query taking so long (only) on first execution?

texts:
I have an AWS lambda function deployed with a spatialite table.  When queried, a point lookup is done on the table to find the locations in the table containing the point.

The first execution of the function I typically see the query execute taking 3.3 seconds.

point_in_polygon_query = 'SELECT ssc_code16, ssc_name16, ste_name16 FROM suburbs WHERE within(PointFromText(\'POINT({lon} {lat})\'),suburbs.GEOMETRY)'
cur.execute(point_in_polygon_query.format(lat=str(lat), lon=str(lon)))

And the fetchall taking 6.6 seconds.
rows = cur.fetchall()

(Determined by using AWS X-Ray subsegments)

However any subsequent runs take about 0.11 seconds each for the query execution and fetchall.

Why is the first execution of my function taking so long (compared to subsequent executions)?

Note: Subsequent query optimisation has significantly improved response times but I still want to understand the why.


texts: How to upload TIGER shapefiles into MongoDB?

texts:
I have the school district shapefiles from TIGER for Michigan, from:
ftp://ftp2.census.gov/geo/tiger/TIGER2012/UNSD/

I'd like to use MongoDBs geoquery features to interrogate this data.

How do I upload the .sh, .sh.xml, .dbf, .prj files into MongoDB so that I can execute these queries?


texts: Why modes are represented with numbers, that are powers of 2.?

texts:
In TraverseModeSet class of  org.opentripplanner.routing.core package, Modes are represented using integer numbers , those are powers of 2. Why?

private static final int MODE_BICYCLE = 1;

private static final int MODE_WALK = 2;

private static final int MODE_CAR = 4;

private static final int MODE_BUS = 16;

private static final int MODE_TRAM = 32;

private static final int MODE_SUBWAY = 64;

private static final int MODE_RAIL = 128;

private static final int MODE_FERRY = 256;

private static final int MODE_CABLE_CAR = 512;

private static final int MODE_GONDOLA = 1024;

private static final int MODE_FUNICULAR = 2048;

private static final int MODE_TRAINISH = MODE_TRAM | MODE_RAIL | MODE_SUBWAY | MODE_FUNICULAR | MODE_GONDOLA;

private static final int MODE_BUSISH = MODE_CABLE_CAR | MODE_BUS;

private static final int MODE_TRANSIT = MODE_TRAINISH | MODE_BUSISH | MODE_FERRY;

private static final int MODE_ALL = MODE_TRANSIT | MODE_WALK | MODE_BICYCLE;

private int modes = 0;


texts:
I would like to do a very simple task: (1) Add a field to a layer and (2) add some values to the field.

I have ADDED A FIELD TO A LAYER by the following code and created three new fields: Dev_site, Buffer and Unit.

Now I want to add values to the field using python. I searched online and found the field calculator tool:

processing.runalg('qgis:fieldcalculator', input_layer, field_name, field_type, field_length, field_precision, new_field, formula, output_layer).

However, in my case, the field calculator is a little bit advanced. I just want to add some very simple string and number to the field which doesn't need any formula.

Is there an alternative (simpler) way to add values to the field using PyQGIS?


texts: Which projection when dealing with least-cost path estimation across the Mediterranean area?

texts: I am planning to perform a least-cost path analysis in ArcGIS 10.1 (advanced license) by means of Path Distance tool. I am taking into account a large area, corresponding to the Mediterranean basin. I am planning to use wind direction as an horizontal factor in the estimation of anisotropic least-cost paths representing navigation routes across that basin. My question: is there any specific type of projection that should be preferentially used when dealing with such a large area and given the type of analysis involved?


texts: Enforce label text within extent of an atlas composer

texts:
I'm not sure if this can be done, but I am using QGIS (2.8) and the wonderful Atlas composer, and what I would like to be able to do, is stop labels from being 'cropped' by the atlas generator (Example below).

I know this goes against the purpose of an atlas creation ( in that there should be overlap so you can navigate from one map to the next), but is there a way to stop this from happening?


texts: How to hide coordinates points in mapinfo?

texts:
I have created labels ,legends scale and North-direction pointer but I now need to conditionally hide or show the points using button in mapbasic program but I don't know the syntax.

  Where can I find syntax examples to
  interact with map window?

I want to change current view in mapinfo window and plot map using values from same table but using certain condition.Any Example.
My Code:

Open Table "geoLocation.TAB" Interactive
Map From geoLocation
select * from geoLocation where vendor = "V1" into selection

Now I want to redraw this selection in same map window but when i do :

Map From Selection

It draws in new window.How to avoid this and is my way correct way ?


texts: Change Feature Type Details in Geoserver

texts:
http://docs.geoserver.org/stable/en/user/webadmin/data/layers.html, Section
"Feature Type Details (Vector)" lists a Layers Features including its Feature Type.

How can I edit these?


texts: Selecting points from intersecting an arbitrary (non-database) polygon with PostGIS

texts:
Determining a polygon as follows, and the points to check against, in ruby

@poly = GEOFACTORY.polygon(@line_string1)
@locations = Location.all
@valid_locations = Location.order("id").where("st_intersects(ST_GeomFromText(@poly), @locations.lonlat)").all

@poly returns  POLYGON ((9.0395777869 45.6141806128, 9.0395777869 45.3141806128, 9.3395777869 45.3141806128, 9.3395777869 45.6141806128, 9.0395777869 45.6141806128))
 and Location.lonlat is defined with :srid => 3785.  Similar queries that check against database column data run properly, however this query states

PG::UndefinedColumn: ERROR:  column "poly" does not exist

Is it possible to query postGIS with arbitrary data?  If so, where is this query off track?


texts:
I am new to both android and QGIS, but have 20 yrs experience with MapInfo and ArcGIS.
I have installed QGIS on a nexus 7 (4.2),  but keep getting the error :

  No qgis data provider plugins found in /data/app_lib/lib No vector
  layers can be loaded.


texts: Find lines not related to any other lines

texts: I have a pipes layer. I need query all dangling pipes, it isn't relationship any other pipes. I have attach image. Is there any method to do this?


texts: How can i use alias of field insteed of field name in arcgis silverlight API?

texts:
I developed a webgis using ArcGIS silverlight API . I'm using ArcGIS 10 and silverlight API 3.1 .I want to show alias name insteed of field name of a feature class. How can i use alias names in my webgis?

Thanks for advance


texts: Trying to create pie charts from points in QGIS show up empty circles

texts:
I'm using QGIS 1.8.

I have a csv file with a list of nature reserves and their coordinates, and there are 10 columns such as 'birdwatching', 'scenery', etc with either 'yes' or 'no' under each for every nature reserve.

What I need to do is display all this as best I can on one map, and since most of the reserves have more than one attraction I figured the best way to do that is with the pie chart function.

When I add in all the attribute fields ('birdwatching', 'scenery' etc) all that I get is empty circles.

I've been playing aroud with settings but can't figure out how to make it work.


texts: Is there a good free Flat style map tile available?

texts: I'm building a website in the style of Flat Design, and hoping to find a map tile that complements the design and can be used with Leafletjs.


texts: Specifying alternate symbology for print vs screen

texts:
Is there a way to specify an alternate symbology for print output in ArcMap (as opposed to what is currently on-screen? I am thinking along the lines of the functionality offered by CSS.

For instance, the on-screen representation of a line might be dashed, but in print it might need to be a double line with a thicker width.

The core issue is that we have users needing to view production layers with overlaid development layers on screen - so a non-standard symbology is helpful here. When design is complete, only the final design layers are printed, but should be represented with the "production" symbology.

Currently using Arc 10.0, but open to solutions requiring other versions and/or scripting.


texts: PostGIS spatial view getting points within distance interval from line

texts:
I have a point table (addresses) and and line table (coastline) in PostGIS. I would like to build a spatial view selecting the points that's has a minimum distance of 600 meters from the line and a maximum distance of 700 meters from the line. I have tried this:

CREATE OR REPLACE VIEW vw_pointsinterval_600_700z WITH (security_barrier=true) AS
 SELECT DISTINCT a.gid,
    a.geom
   FROM coastline c
     LEFT JOIN address a ON st_dwithin(a.geom, n.geom, 700::double precision) AND not st_dwithin(a.geom, n.geom, 600::double precision);

The spatial view is visualized in QGIS added two test buffers. One of 600 and one 700 meters from the coastline. This does not show the correct result:

The coastline is the red line, the purple color is a test buffer of 600 meters from the coastline. The light green is the 700 meter buffer and the green dots are addresses. Some random places the view returns addresses that's closer than 600 meters to the coastline? I would expect all the green points to be on the green polygon test buffer.


texts: Error creating DEM with Nearest Neighbour tool in QGIS?

texts:
I have an issue when I am creating a DEM with the tool natural neighbour on QGIS.

When I extract my data directly on Qgis I have no issue creating a DEM. My points are coming from different sources and I need to join them all together.

When I use the tool Union and then create my DEM with those points I have the following error message:

I don't know what it means and I tried to save the points directly on C: with no space.

It seems to be a minor issue but can't find a solution


texts: Changing map scale in ArcMap after almost finished digitizing

texts:
I'm just about to design a map that represents different areas in my immediate living environment. I Unfortunately made a big mistake at the beginning of my project. The current scale is given by ArcMap with 1:30,000. 1: 3,000 would be the correct scale.

Is it possible to transform the already almost finalized digitization in any way to the correct scale, so that the distances as well as the areas in my map correspond to the reality?


texts: Flow Distance Not Calculating for Entire Raster

texts:
I am trying to calculate the height above the nearest stream for an area. The problem is, it only does so for a portion ( the colored area below) of my area of interest. In the pic below, the black and white raster is my DEM that all of the work is based off of. I want my results of the Flow Distance tool to cover that entire DEM. It missed the lower smallest stream entirely, and the lower portion of another.

In Environments > Processing extent and Raster Analysis mask, I have both set to the DEM. How come it won't calculate a value (height above nearest stream) for the entire area?

My inputs to Flow Distance are a raster of the 3 streams shown below, surface raster is the DEM (filled), and I have tried both with and without a directional DEM, matching it correctly to how it was made (D8).


texts: Calculate buffer / multi-ring buffer of geometry with GDAL / gdalUtils

texts:
I wonder if there is any way to calculate a buffer on a spatial data set with basic command line GDAL or the R interface gdalUtils.

To specify the Q: I have a lot of points around the globe. For these points I have to calculate (multiring) buffers. Since this should take place in an automated framework implemented in a mix of R and bash, I wonder if there is any way to use the power of GDAL/OGR?

I found this geometry class during my research on this problem. However, I am not sure if this can be part of the solution.

For test purposes, let's assume the point is lon/lat (EPSG:4326): 0, 0.


texts: Transferring Interpolated Results to Line Features

texts:
I am using ArcDesktop 10.3.  I have a road network feature dataset and a point city feature class. Using Network Analyst and Closest Facility along with several incident points, I solved several routes for distances to the incidents. I then joined the distance from facility to the incident to the incident.  I then used IDW to interpolate travel distance to all areas on my map.

I like how it looks, but I had an idea to transfer the color of the travel time (the light to dark blue in the map) not on the country, but rather on the road (line) itself.  So the road as it leaves one of the "facility" cities, it will start as light and become darker the further along in distance it goes.

I thought I would join the raster value I created in the IDW to the line, but I can't seem to do this.  How would one go about accomplishing this task?


texts: FME: Set attributes based on a lookup table

texts:
I have an XML Reader and I need to write (both inserts and updates) to an SQL database.  The attribute data  in the XML is in a string format, however the values in SQL are stored as IDs.  There is a reference table in SQL that stores a lookup between the field name (ComboType in the screen shot below), string value (Code), and the stored ID (ID).

So for this example, I need to insert/update 5294 when the "Digger" value is given for the SW_Equipment_Needed tag.  I know that the there is some mapping required between XML and SQL.

What is the best way to retrieve the stored ID value for each of my feature attributes?  There could be any number of attributes per feature depending on what is in the source XML.  There are over 6,000 of these lookup values.  My first instinct is to use Python and the pyodbc module but there could be a cleaner way.


texts:
I have a lot of coordinates I want to convert from Sweref99tm to WGS84 for display on google maps. I'm using gdal installed on my laptop:

$ echo "6782384 599140" | gdaltransform -s_srs EPSG:3006 -t_srs EPSG:4326
64.1363442998149 3.53655177487795 0

But according to several sites online where I input the same coordinates, I get 61.163396, 16.842654 returned instead.

Anyone have an idea what I'm doing wrong?


texts: arcpy clean up of variables in cursors vs. loops

texts:
I'm new to Python & ArcPy. I'm used to cleaning up temporary variables used in loops after completing the loop to keep my workspace clean. I think I read that variables created in cursors are temporary and deleted once the cursor completes? In the code sample below, what (if any) should my clean up steps be?

counts = {}
with arcpy.da.SearchCursor (outFeatureClass.shp, "Zone_ID") as cursor:
    for row in cursor:
        if not row[0] in counts.keys():
            counts[row[0]] = 1
        else:
            counts[row[0]] += 1
# del row?? cursor??

with arcpy.da.UpdateCursor(outFeatureClass.shp, ["Zone_ID", "N_1"]) as cursor:
    for row in cursor:
        zone = row[0]
        if zone in counts:
            row[1] = counts[zone]
            cursor.updateRow(row)
# del row?? zone?? cursor??

with arcpy.da.UpdateCursor(outFeatureClass4.shp, ["Zone_ID", "N_1"]) as cursor:
    for row in cursor:
        zone = row[0]
        if zone in counts:
            row[1] = counts[zone]
            cursor.updateRow(row)
# del row?? zone?? cursor??


texts: I'd like to know if there are some maps giving the geographical coverage of the media publication, for example, a map giving the counties or even the cities where local newspapers are sold, and a map showing the coverage of the US radios?


texts:
Context

I have a database with spatial data such as areas and points. Some of the areas in my db are coming from OpenStreetMaps, representing some regions. (Which are kinda small at earth's scale, but wide enough to be problematic)

I'm using PostGIS extension to keep them in my db, as GEOMETRY objects (with 4236' id). They're coming with their latitude / longitude values, so they're finally just a succession of points, geolocated. So, my current objective is to project these objects on a map (such as GoogleMap), there come the problems...

Google Map Overlay

I use a custom "overlay" object. Creating my overlay with the svg's envelope (got from ST_Envelope PostGis function). Then, in this "bounding rect", I will draw my svg.

This bounding rectangle is pasted as a div on the map, with the google functions "fromLatLngToDivPixel(bounds)". The div comes then exactly where it has to be drawn.

Drawing the SVG

The bounding rect is the "svg container" on my map. So, inside of it, I come and paste my svg shape. First, I build a "svg object" with attributes that will make it "fit" the div and make the contents (path) visible, such a the viewbox.

For a geo-object coming from 50° North and 2° East and smaller than 1° both sides, then my viewbox looks like : "viewBox='2.00 -50.00 1.00 1.00'".

After that, I have to consider that the container (div) is not at the same scale as the svg' object, so, my current "solution" which is not working is to "scale" my svg to fit all space in the container div : "preserveAspectRatio='none'" AND "style='width:100%;height:100%;'".

At this point, I was about to say "Hell Yeah! I did it!" but...

Problem

The svg' points are placed "as if they were directly placed on a lat/lon system" but what I'm actually doing is draw it as a svg on a "2d linear system". Meaning that the distance between 2° east and 3° east is the same as the distance between 3° east and 4° east WHEREVER I'M LOCATED NORTH OR SOUTH! This assumption, which is completely false, brings me to that kind of diplay problem :

At some place of the svg, the boundary I'm drawing perfectly fits the map's boundary. BUT at other places, a small shift is clearly visible on the map, due to this "scale" problem.

Note that, at the very bottom of the map, my boundary is correctly placed! but the more you go north, the more it's shifted...

Concrete example on jsFiddle : http://jsfiddle.net/jLw2umt6/

How to solve?

Couldn't find a decent solution yet, already thought about these :

Change from svg to "image" such as png and so on with help of a function like the "fromLatLngToDivPixel" used to place the bounding rectangle...
Divide the shape into many many smaller shapes, placed into several smaller bounding rect's (overlays), but performance issues... :/
Find a function somewhere, somehow, it's probably already available somewhere, that will, instead of returning the svg as it is now, will return "projected points" based on the svg lng/lat location...
Find another way (maybe directly from postgis object or from another tool) to project a shape, area on google map?
Find a way to make earth cubic, or completely flat instead of spherical as it is now.


texts: Seeking package/function in R that can help tag GPS coordinates of fauna sightings to nearest sampling point?

texts:
I have 1 dataset of animals sighted walking a survey, and I would like to tag them to the respective sampling points along a transect.

Any idea what function/package in R can do that? Maybe something nearest neighbour related?


texts: How to get around the 10-character limit when converting from kml to shp with ogr2ogr

texts:
I'm trying to create custom .shp files by drawing polygons on a new layer in Google Maps Engine, exporting that .kml file, and then converting that file to a .shp file with ogr2ogr.

If I run a command like this:

ogr2ogr -f 'ESRI Shapefile' foo.shp foo.kml

I get the following printed to stdout:

Warning 6: Normalized/laundered field name: 'Description' to 'Descriptio'.

This causes a problem downstream because when I load the .shp file into a r-tree and then query it to perform reverse geocoding, I get back a field "Descriptio" instead of "Description", which just looks odd. Moreover, this happens despite the fact that when I open the .kml file in vim, I don't see any "description" tag anywhere.

So my question is: is there any way around this? Can I alter the KML file in any way to prevent ogr2ogr from complaining that it needs to normalize/launder the field Description?

Please let me know if this question is written poorly. It's my first time on GIS Stack Exchange and I'm not at all trying to be difficult/annoying/insufferable.


texts:
could I change the icon of my legend in the mapcomposer?
I want other icons for my layers in the legend then in the project.
I use the Version 3.4.8 LTR.