texts:
I'm looking for the best ways to 'group' a budge of block polygons in a new polygon that should be the town.

My block layer is something like this:

Or zoomed:

And I'm looking to get something like this polygon (the left one):

Could somebody tell me some indications of how to get something like this?
I'm not looking for the best accuracy, but something close to town or city polygons.

I thought in a convex hull, but the problem is that the block layer are a lot of polygons, instead of an irregular one.

I'm sure that PostGIS has some process that approximates to my goal, but I can't find it.


texts: How to delete a file after processing an algorithm in QGIS 2.18?

texts:
I am trying to delete an intermediate file ("Temp_merged.shp") after running the singleparts to multiparts algorithm in QGIS 2.18.

   import glob, os, processing
   path_to_shp = "C:/Temp/Tracks/"
   os.chdir(path_to_shp)

   for fname in glob.glob("Temp_merged.shp"):
        merger = "Temp_merged.shp"
        processing.runalg("qgis:singlepartstomultipart",merger,"Unique","Multipart.shp", progress=None)

   del merger
   os.remove("Temp_merged.shp")

However, I am getting the following error:

  [Error 32] The process cannot access the file because it is being used by another process: 'Temp_merged.shp' See log for more details

How can I end the processing.runalg to delete the intermediate shapefile?


texts: How to Choose a Projected Coordinate System for Nepal?

texts:
I am working on a project mapping the tiger habitats in Nepal (predominantly along the southern border of the country) and I am not sure which coordinate system / projection I should be working in.

Any suggestions?


texts: Generating exportable report from web map using ArcGIS API for JavaScript?

texts:
I would like to build a web map using ArcGIS JavaScript API containing 1 layer of polygons & 1+ layers of points (or lines) contained within said polygons. I want this to be enabled with some sort of exportable "report"-generation functionality, which can do the following:

In exportable tabular format, show the points/lines contained within each polygon and their attributes, i.e. the result of a one-to-many spatial join between the polygons and points.
Optional for the user, a table like in #1 but aggregated by some summary field of the points/lines (specified by the user via drop-downs menus or whatever), like the tabular part of a Dissolve operation in desktop ArcGIS, on some specific dissolve field. So for each polygon, they could read the total length of the lines it contains, etc.

The report table doesn't have to show up in the browser (though that would be quite desirable), but does have to be exportable to some common format like CSV or preferably Excel.

Does anyone know the best/quickest way to accomplish this?

I know just enough about JavaScript/HTML/the API to have ideas to pursue but not enough to really know which are the most viable.

Some ideas:

Reading table attribute values directly via a queryTask & printing them out (à la this ESRI sample), but in a more complex tabular format, e.g. dynamically filling a nice-looking HTML table with the attribute values, then exporting the table as CSV (which cursory Googling suggests might be possible in one way or another).
Making use of a dgrid at least for display, but somehow linked to the above table manipulations. I don't know much about dgrid. Then export this, which appears to be possible with the Exporter plug-in for Dojox EnhancedGrid.

Has anyone implemented this kind of reporting functionality via an ArcGIS JavaScript API web map?

Is there any easier way to do this or other possibilities?

Any warnings about the viability of the above?

I just want to scope out a smart method before striking out blindly.


texts:
How to remove some icons (not all) on the Layers Panel dockwidget using python? Icons on the Layers Panel are :

I have written the code :

from PyQt4.QtGui import QDockWidget
for icon in iface.mainWindow().findChild(QDockWidget,'Layers').actions():
print icon.objectName()

with this code, nothing happen. So, I can not remove certain icon like in Name of children of toolbar in QGIS python and how to disable?


texts: Lowrance sl2 format... is it GIS data and/or can it be made into GIS data

texts:
I have been asked to investigate if Lowrance SL2 (side-scan sonar) format can be put into maps for reports.

As I understand it the SL2 format comprises 3 "bands" left-down-right and is more of a video than raster. If the sonar log file has been captured previously can it be opened and converted to a raster and then viewed - prefer with georeference?

The current thoughts are to screen-shot the log using the sonar log viewer and then georeference but having seen the contents of these files I don't like the chances of getting similarity enough to georeference.


texts: Export Table to Excel works in ModelBuilder but not in MXD?

texts:
I have created model that take user input from a picklist and then outputs the selected data into an Excel File. When I run the model inside of ModelBuilder it outputs an Excel file that is exactly as it should be. Then, when I save the model and then close ModelBuilder and try to use it in an MXD, the model runs successfully and then give me a path to the output (the very same path that ModelBuilder gives) but when I follow the path there is no spreadsheet.

Has anyone experienced this before?


texts: What is the latest version of R that supported gpclib?

texts: Does anyone know what the latest version of R that supported the gpclib package? I was using the latest version but when I tried to install gpclib, it said that it was unable due to my version of R. The gpclib page gives the minimum version and requirements. Do I have to downgrade my version of R?


texts:
My objective is to convert a google earth engine featureCollection into a geopandas geodataframe.

Earth engine can handle geoJSON whereas geopandas prefers shapely.

The geometry (geoJSON) of the first feature in my featureCollection is:

https://pastebin.com/FSWFpBu4

as you can see, there is a geodesic (boolean) parameter in the geoJSON.

I use shapely.geometry.shape() to convert the geoJSON to a shapely geometry that geopandas can handle. However I am not sure if the geodesic parameter is captured correctly. The shapely documentation doesn't mention this. [EDIT: the main page does] How can I make sure the geodesic parameter is handled correctly?

Here is my full function:

def fcToGdf(fc):
    """converts a featurecollection to a (geo)Pandas GeoDataFrame or a pandas DataFrame. WARNING: Geometries are non-geodesic. Geodesic not yet supported. work in progress

    Args:
        fc (ee.FeatureCollection) : the earth engine feature collection to convert. Size is limited to memory (geopandas limitation)
        crs (dictionary, optional) : the coordinate reference system in geopandas format. Defaults to {'init' :'epsg:4326'}

    Returns:
        gdf (geoPandas.GeoDataFrame or pandas.DataFrame) : the corresponding (geo)dataframe.

    """
    crs = {'init' :'epsg:4326'}

    features = fc.getInfo()['features']
    dictarr = []

    for f in features:
        attr = f['properties']
        attr['geometry'] = f['geometry']
        dictarr.append(attr)

    gdf = gpd.GeoDataFrame(dictarr)
    gdf['geometry2'] = map(lambda s: shapely.geometry.shape(s), gdf.geometry)
    gdf.crs = crs
    return gdf


texts: How to convert dbf to raster the short way?

texts:
I have got a dbf file and I want make it to raster file but the short way,

Look at my long way,

First i add my dbf file
and i create tin
http://i1211.photobucket.com/albums/cc421/mtnkou/2.jpg
at the and my raster file (tin to raster)
http://i1211.photobucket.com/albums/cc421/mtnkou/3.jpg

How can I do  dbf to raster with short way,

i don't want create tin file,

Just dbf to raster,

How can i do ?


texts: Compute the distance of some nodes using a road layer which are not exactly on the road

texts:
I have the road network of a region (downloaded from open street map as a layer) and a set of nodes (selected using google earth, saved as a KML file and finally converted as a layer).

I need to calculate the shortest distance between the nodes using the road (not the Euclidean distance!), remember that some of the nodes are maybe far from the road and not exactly located on. The output should be a matrix distance showing the distance for each pair of nodes.


texts: Converting DEM text file to tiff file in QGIS?

texts:
I am trying to convert a digital elevation model in txt file format to a raster file (tiff) with qgis.

The file format is *.dat. The coordinate system is Gauss Krüger, Bessel, 5. meridian.

First I imported the txt file with the appropriate coordinate system. Now I want to convert it to a raster file. However the viewed data are black and I dont see any heights. What did I do wrong? Follow the link to view the screenshot of the imported dem in qgis:
Screenshot qgis

When zooming in it looks like this:
Screenshot zoomed in

This is an example of the data:

LATITUDE    LONGITUDE   HEIGHT
5390.0125   5659.9875   258.6
5390.0125   5659.9625   252.6
5390.0125   5659.9375   243.9
5390.0125   5659.9125   237.1
5390.0125   5659.8875   234.5
5390.0125   5659.8625   231.2
5390.0125   5659.8375   233.7


texts: QgsExpression using complex sql and python subsitution

texts:
I have a sql statement i am running in postgres fine and would like to transfer it to a QgsExpression mainly because I have loaded the layer in Qgis.
My sql is:
 select column a,column b from
table1  GROUP BY column a,column b ORDER BY column b DESC limit 20.  I do not know how to do so. Can I use python substitution when creating an expression.

expr = QgsExpression( "%s = 99" )% column a


texts: Batch merge large numbers of.asc grid files

texts:
I have >2.5k .asc grid files which I need to merge into a single overarching dataset.  I would like to avoid having to convert them into esri grids first so was looking to use QGIS (v.2.0).

I am aware of the command line length limitation with the QGIS merge tool so was looking to use the OSGeo4W shell to automate the process.  Can someone help me with some code to cycle through the files and merge them?

I was looking at this thread but the suggested solution doesn't work - it results in an "index out of range error"

The files are all of the same format, in a single folder and have the name convention SDXXX.asc - where X is a single digit.


texts: QGIS 3 Installation fails - Windows 10 - Main Window not responding

texts:
I am trying to install QGIS on a win 10 Desktop PC, on which QGIS 3 was running until Win 10 has been reseted and compeletly new installed. After the installation, the QGIS splash screen shows up (QGIS is ready) but the main window remains empty (-not responding-, as reported by some other users as well).
I have tried allmost all versions 3.10, 3.14, 64/32Bit and Installation with OSgeo4w - always the same result: QGIS not main window is not responding. I also have removed local Appdata before every installation, also all hardware drivers seem to be up-to-date.
Any Ideas or Suggestions? Are there any logfiles available, which might help to find the Problem? Are there any additional tools/ programs needed?


texts: Creating 700 m x 700m grid in ArcGIS for Desktop?

texts:
I'd like to know how to create 700m x 700m grid array in ArcGIS.

I am a beginner at GIS.

I'd like to develop the insect distribution map.


texts: Determine the geotransformation to convert a netCDF to geotiff

texts:
I have a netCDF file I need to convert to a geotiff.  There are three subdatasets in the netCDF file: snow, lat and lon.  I am accessing these files like this:

import gdal
import netCDF4 as nc

ncfile = nc.Dataset("C:/path_to_file/file.nc", "r")

snow = ncfile.variables["snc"]
lat = ncfile.variables["lat"]
lon = ncfile.variables["lon"]

The metadata says the latitude and longitude are in degrees, so I want to project to WGS.

The minimum and maximum latitude and longitudes found by:

xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]

returns:

(0.0, -87.863801343710804, 357.1875, 87.863801343710804)

since my longitude goes from 0 to 360, and it should go from -180 to 180 I am adjusting it like so:

lon = ((lon+180) % 360) - 180

I now determine my geotransform:

nx = len(lon)
ny = len(lat)
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)

and try to create my tif:

dst_ds = gdal.GetDriverByName('GTiff').Create('output.tif', ny, nx, 1, gdal.GDT_Float32)

dst_ds.SetGeoTransform(geotransform)    # specify coords
srs = osr.SpatialReference()            # establish encoding
srs.ImportFromEPSG(3857)                # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(snow)   # write r-band to the raster
dst_ds.FlushCache()                     # write to disk
dst_ds = None

I followed this unanswered question to get this far, https://stackoverflow.com/questions/43254024/extracting-specific-netcdf-info-and-converting-to-geotiff-in-python, but my raster won't even output.

I instead get this error:

Traceback (most recent call last):

  File "<ipython-input-10-0384d872ba66>", line 15, in <module>
    dst_ds.GetRasterBand(1).WriteArray(snow)   # write r-band to the raster

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda_64\lib\site-packages\osgeo\gdal.py", line 2332, in WriteArray
    callback_data = callback_data )

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda_64\lib\site-packages\osgeo\gdal_array.py", line 365, in BandWriteArray
    raise ValueError("array larger than output file, or offset off edge")

ValueError: array larger than output file, or offset off edge


texts: Downsampling then upsampling raster using QGIS?

texts:
I'm trying to downsample a raster then upsample it again.

I am using the gdalwarp function in QGIS.

When downsampling the raster looks OK in the front end, I set the cell size to 10x the original raster.

When upsampling back to the original cell size using the downsampled raster the original data seems to be preserved and I end up back at the original raster. What I'm expecting to happen is that I see blocks of 10x10 cells all with the same value (same as the downsampled raster).

Does QGIS preserve the original cell values?

Is there a way to perform this operation?


texts: WFS, ArcGIS Server 10.2 and WFS maxFeatures?

texts:
I'm experiencing intermittent 'out of server memory' errors on
our public ArcGIS Server 10.2 server.

From what I can glean from the OGC documentation the server side
maxfeatures parameter was only available in the WFS 1.1.0 spec but
not 1.0.0

I suspect that our server is temporarily being overwhelmed by WFS 1.0.0 requests.

IE some feature types have in excess of 90,000 records

We need WFS capability on our services for integration with another system. The other
system is being coded to perform WFS 1.1.0 requests so the maxfeatures will apply.

Is there any way I can 'turn off' WFS 1.0.0 support in ArcGIS Server 10.2 without turning
1.1.0 off?

The only thing ESRI support and I came up with was to create a custom getcapabilities
xml version 1.0.0 file and delete all feature types in that minus the one with least records...

Regards,
Richard


texts:
I have a largish raster file (100 MB in size) with resolution of 0.05 degree. This raster file has temperature data for the entire world, but I only need data for the U.S. I have one such file for each day of the year from 2000 to 2015.

Is there a way to save time while reading this file using rasterio perhaps by specifying the extent of U.S?


texts: AttributeError: Result object has no attribute save (when trying to save output in separate directory)?

texts:
I am trying to convert more than 100 ascii files to raster using python code. All the lines (except Last two) in my code presented below is successfully converting Ascii to tif format in the working directory. However, I want to save the new TIF files in a separate directory. Therefore, the last two lines in my code, I am trying to save in a new directory. But, I am getting an error:

"Traceback (most recent call last):
  File "C:\Subhasis\Project-Biological_monitoring\R_SAGA\Raster_data\ASCII_Raster.py", line 22, in <module>
    Ascii_raster.save(outname)
AttributeError: 'Result' object has no attribute 'save'"

Code:

import arcpy, os
from arcpy import env
from arcpy.sa import *

#To overwrite output
arcpy.env.overwriteOutput = True

#set work environment
env.workspace = "C:/Subhasis/Project-Biological_monitoring/R_SAGA/Raster_data"

outWorkspace="C:/Subhasis/Project-Biological_monitoring/R_SAGA/Raster_data/Final"

dataType="INTEGER"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

for file in arcpy.ListFiles("*.asc"):
    outputraster = file.strip(".asc")+ ".tif"
    Ascii_raster=arcpy.ASCIIToRaster_conversion(file,outputraster,dataType)
    outname=os.path.join(outWorkspace,outputraster)
    Ascii_raster.save(outname)


texts:
How can I install new version of QGIS (2.0)?

When I try with the indication of the official web site of QGIS, terminal write me that the last version of QGIS (1.7.5 wrowclaw) has several plugins that not permise it.

I am a new user of ubuntu 13.10.


texts: This should be an easy task but I have struggled for 2 days trying to do it. I have a Qfield project with a single point layer (set to offline edit). I have a simple basemap (OpenstreetMap). I collect data points on the handheld and all looks good. I copy the whole folder containing the project back in to the import folder for QGis sych. When I import the project the points are missing. What am I doing wrong?


texts: How to calculate the area of a polygon given its coordinates?


texts:
I have a problem, but I didn't find usable solution. I created a line shape as a reference/guideline and made grid with grass v.mkgrid and I should correctly fit the grid onto the line.

It isn't possible to fitted precision with manual method, and QGIS Geometry snapping didn't work correctly.I approached the line with move and rotation, the distance within 5 meter

If I used geometry snapper, the tool modified shape and size of polygon, connected to line.
Other situation, if line was within the grid, the tool worked as clip.

Somebody knows such QGIS tool/plugin  what solving this problem?

There is nothing logic between layers, line attributes contains just one id, grid contains mkgrid generated data

I tried snap geometry to layer, the results was wrong, the tool stretched the cells to the line, modified  first column's size and because the grid wasn't parallel to line formed the cells into trapeze


texts: How to create a date field in QGIS shapefile

texts:
When I create a new shapefile in QGIS, I'm presented with the "New Vector Layer" window, where I specify the Shapefile type (Point/Line/Polygon) and am then able to create attribute fields.

I'm only able to select the field type as Text/Integer/Decimal - how do I go about creating a date field type?


texts: Styling GeoServer vector tiles using *.sld files (PostGIS, QGIS, Leaflet)

texts:
I've been styling vector layers published to GeoServer using *.sld files saved from QGIS as explained in this link.

Now, I'd like to work with vector tiles. While I am able to code the JSON within the vectorTileLayerStyles parameters, I am wondering if there is any way to use the styles set from the SLD files instead.


texts: Alternatives to ArcGIS World Geocoder for searching coordinates in AGOL web app

texts:
I'm building a web app in ArcGIS Online, and don't really care for the way the ArcGIS World Geocoder (https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer) handles coordinates for the search widget.  Are there any alternatives that might be a bit more user friendly for end users of the app?

Currently, I'm using the coordinate conversion widget as a workaround, since it has "zoom to" and "add point" functionality, but would like to know if there are decent geocoder alternatives for searching coordinates from the search widget (both DD and DMS, if possible).


texts: Authenticating ArcGIS Server Feature Service in QGIS?

texts:
I am trying to add a protected Arc layer to a QGIS project and have had no success.

What I have tried so far:

1) Generate token using the token service url...

2) Create a QGIS authentication configuration...

3) Enter the service url...
https://geodata-services.drillinginfo.com/arcgis/rest/services/Landtrac/DI_Landtrac_Units/FeatureServer

When attempting to add the layer I get...

Failed to retrieve service capabilities:
Error 499: Token Required

However, if I append the token to the url with: url?token=mytoken, I am able to see the layer in the list...

But I am unable to add the layer to the map...
I receive error:

2019-09-26T12:52:55     CRITICAL    Layer is not valid : The layer authcfg='10ph3f0' bbox='-1.13707e+07,3.75972e+06,-1.13549e+07,3.76784e+06' crs='EPSG:3857' filter='' url='https://geodata-services.drillinginfo.com/arcgis/rest/services/Landtrac/DI_Landtrac_Units/FeatureServer?token=mytoken/0' table="" sql= is not a valid layer and can not be added to the map. Reason: AFSProvider getLayerInfo failed


texts:
I've been reading through the help and readme.md and playing with the configuration options for the CMV LayerControl and I'm pretty confused by the possibilities. It seems to me there are at least 2 if not 3 places where configurations can be set for this control.

The first choice seems to be where you configure the widget in viewer.js in the widgets.layerControl.options. In the "options" object you can add "Control Options" configuration attributes such as "exclude", "noLegend", "noZoom", "noTransparency", "swipe", "swipeScope", "expanded", "sublayers", "metadataUrl", "allSublayerToggles". Setting it here would affect all service "layers".

The next place where you can do some configuration is the viewer.js in the operationalLayers[].layerControlLayerInfos object. It looks like here you can also add some "Control Options" at the service "layer" level. This would affect each service "layer" independently.

The final place, and this is what I'm most confused about, it appears that you can create your own LayerInfos array of LayerInfo objects to customize the Layer Control for each GIS "layer" (ie. sub-layer in a Dynamic Map Server Layer)?

Maybe I'm just wishing this to be the case.

Can someone confirm that?

Where would you use the LayerControl Class?

I guess that would be if you wanted to add the LayerControl widget to the DOM in a some other location?

Can you do anything with sub-layers in a Dynamic Map Service Layer?


texts: Map Info: Thematic Mapping Symbols part printing

texts:
I have produced a series of thematic maps for graduated housing figures, using a symbol which I customized, in three different tables displayed in a layout.

This looks fine in the map window, although the symbols are sized much larger than in the layout and do overlap, in order to be a reasonably large size in the layout window.

All looks fine in the layout window, no overlapping. However, when printed, a couple of the symbols have parts missing. One has the top half and the other the right hand half. The symbols for these tables are small dots and don't overlap.

Any ideas please?


texts: How could I insert a 3-D LineString into PostGIS table?

texts:
Using:

CREATE TABLE trajectories(
gid SERIAL PRIMARY KEY,
vehID VARCHAR(10),
state INTEGER,
geom geometry(MULTILINESTRING)
)

I get an error, for instance, after inserting:

INSERT INTO trajectories(vehid,state,geom)
VALUES(001140,2,ST_GeomFromText('MULTILINESTRING((1 2 3,4 5 6))'))

it gives:

  Geometry has Z dimension but column does not.

How could I circumvent this error?


texts: Publishing map to ArcGIS Server with time enabled data and query joining transposed data with shapefile

texts: I am trying to figure out if I can publish my map service to ArcGIS Server that has time enabled data and uses a query table.... The query table joins my transposed time data and my shapefile. When I tried publishing it, I received the error "00004 layer's data source is not supported" and have not been able to figure it out....I have my shapefile and the transposed data in a file geodatabase....I am hoping there is some quick fix for this, but if this is not possible, is there some other option?


texts: Google Earth Engine extract data from Sentinel-5P

texts:
I would like to extract Near Real-Time Nitrogen Dioxide from Sentinel-5P for all countries. The output (in csv format) however does not show the Nitrogen Dioxide information and date (only country ID is displayed). The world shapefile can be found from: https://code.earthengine.google.com/?asset=users/tronganhtrinh/worldsp
// Define which year to work with
var yr = 2020;

// Range for the Landsat images
var Range = ee.DateRange(yr + '-05-28', yr + '-05-31');

var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
  .select('NO2_column_number_density')
  .filterBounds(worldsp)
  .filterDate(Range);

function tabulate(i){
  return worldsp.map(function (f){
    var r = i.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: f.geometry(),
      scale: 500,
      bestEffort:true,
      maxPixels:1e9});
    return ee.Feature(null, {
      name: f.get('OBJECTID'),
      index: i.get('system:time_start')});
  });}

var resultstable = collection.map(tabulate).flatten();

Export.table.toDrive(resultstable,'results','csv');


texts: List all Unavailable layers in QGIS (3.x) project with PyQGIS

texts:
I am looking for a Python macro to list layers which are unavailable in a QGIS project with a large number of layers.

The following code lists all layers in the current project:

from qgis.core import QgsProject
layers = QgsProject.instance().mapLayers()
print(layers)

I would like to also show the status (available, unavailable) for each of the layers. How can I do this?


texts: Summarize table by group and add logic to keep data

texts:
I need to dissolve a polygon layer by group but, at the same time, I need to keep some of the data in extra columns. This is the diagram:

This is the logic for the new group:

Cod = "Cod" where N1 is the highest
N1 = highest "N1" of group A
N2 = highest "N2" of group A
N3 = get value of N1 where "Cod" is the longest in length characters
N4 = get value of N1 where "Cod" is different from 0 AND 1 AND 1_2

The issue here is that I need to incorporate this into the modeler. Do I have any possibility. By the way, I could use either QGIS 2.XX or 3.XX


texts: ImportError installing GDAL on Ubuntu 14.04 with Anaconda Python 3.5

texts:
I installed the GDAL library on Ubuntu 14.04 through the Anaconda (Python 3.5) distribution. After installation, which seems to work fine, I try importing it and I get this error:

>>>from osgeo import gdal

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/klab/anaconda3/lib/python3.5/site-packages/osgeo/__init__.py", line 21, in <module> _gdal = swig_import_helper()
  File "/home/klab/anaconda3/lib/python3.5/site packages/osgeo/__init__.py", line 17, in swig_import_helper_mod = imp.load_module('_gdal', fp, pathname, description)
  File "/home/klab/anaconda3/lib/python3.5/imp.py", line 242, in load_module return load_dynamic(name, filename, file)
  File "/home/klab/anaconda3/lib/python3.5/imp.py", line 342, in load_dynamic return _load(spec)

ImportError: libcom_err.so.3: cannot open shared object file: No such file or directory

I tried installing through PIP and Conda install commands with the same error. I did go to the directory where this should exist, and it indeed is missing. I figured that if this were a true dependency, it would have installed automatically. Has anybody else had this problem?

I was trying to follow this link, but it wasn't helpful. https://groups.google.com/a/continuum.io/forum/#!topic/anaconda/2-bXTbSiQzg


texts: Selecting points within fishnet grid cell with criteria of at least 3 unique values using ArcGIS Desktop?

texts:
My goal is to create a raster from points with values based on at least 3 unique IDs within a cell.  Cells that had fewer than 3 points with unique IDs will have value NA. See below for clarification. I looked at Selecting set number of points within polygon grid cells using QGIS? hoping to adapt it, but am still stuck.

I created a fishnet with certain cell size that I think may be reasonable. Now, I want to identify the number points within each cell with unique IDs (here, VesselID). Next, I want to assign NAs to my field of interest (here, AdjPounds) for points that do not fit the "minimum 3 unique IDs in the cell criteria". I thought about using the Intersect tool, but don't know how to add the criteria.

Can you help?

This is probably best done in Python, but I'm not familiar with that language yet.

More context: I will need to reiterate this process for several similar files and ultimately conduct raster math to get an average raster. I assume I will need to choose the same cell size when I create the rasters to avoid resampling to obtain the final "average" raster.


texts: How to program this Delete case in QGIS?

texts:
Oke I think I got the hang of cases now, but still I am learning. So I kind of am still making mistakes and I cannot exacly understand why.

Currently I want certain value's in a certain column to be deleted.

I am adding the field "Delete" expression to the column in edit attribute table. So I assume I do not need to add the column name since I am setting the Field calculator in that column:
This is the expression:

DELETE '2013' and '2015' and '2016'
END

However I always get syntax error, I am assume I am programming the expression wrong? Please let me know.


texts: Multiple layers within a single geojson file

texts:
To give some background, I'm working on a project that has quite a few restrictions given the environment it will run in. I need to be able to view and interact with geospatial data, perform geospatial data entry (ie click a point on the map and fill in the attributes), no new software or server can be installed on the client machine, and it cannot be added to an existing web server.

My solution is to use OpenLayers from an html file opened off of the client's file system. This resulted in cross-origin security violations in browsers when attempting to read spatial data files off of the local hard drive from javascript (despite all the files being from the same directory), and modifying the browser security is not an option. My workaround for this is to add a file drag and drop where the user "uploads" the file and it's read into OpenLayers, but I have to make it easy for the user and dropping 10+ files individually isn't practical.

I want to merge each of the 10+ GeoJSON files I use into a single file, but keep them distinct, such that I can load them into OpenLayers as individual layers in order to apply individual styling and not try to merge or generalize unrelated attributes. I can't think of a way to do this cleanly where OpenLayers will still be able to interpret the GeoJSON structure.

Is it practical to wrap the GeoJSON in XML and break it apart in javascript before sending it to OpenLayers? Can I wrap GeoJSON in JSON or would it confuse parsers?


texts:
I have a shapefile:

Street1    Street2    Street3
Field      Clark      Taylor
Field      Jackson    State
Franklin   23rd       25th

And a Dataframe df:

Page    Entry_Date
12      08012018
17      07252018
11      07252018

The only commonalities between the two are the index position of the rows. I'd like to join them together so my result is as below in my shapefile

Street1    Street2    Street3    Page    Entry_Date
Field      Clark      Taylor     12      08012018
Field      Jackson    State      17      07252018
Franklin   23rd       25th       11      07252018

I've tried a couple ways but it doesn't seem that arcpy likes to cooperate with pandas

in_shp = "C:\\Users\\Anthony\\Proj56\\ExportShp.shp"

for idx, row in df.itterrows():
    arcpy.JoinField_management(in_shp, "FID", df, idx, [df.columns.values])

This just gives me an error. Is there any way to incorporate a pandas dataframe into an existing shapefile without completely rewriting everything?


texts: Extract first and last element of path data type in Postgres?

texts:
Our postgres database has tables that have a "path" data type consisting of a list of points that make up a line. Is there any way in postGIS or postgres to extract the first and last coordinate pairs from each entry? The purpose is to determine the cardinal direction of each line using ST_Azimuth, but that function requires the coordinates of the first and last point.

I've tried every way I can think of, but the path data type doesn't act like a string or an array, so it's difficult to pull out items from it.


texts: How to remove 'no data' values in sentinel 1 data? - Sentinel Toolbox

texts:
I hope anyone can help with this issue:
I am using the Sentinel Toolbox to preprocess Sentinel 1 GRDH scenes. Managed to radiometric calibrate, conduct radiometric terrain flattening and geometric terrain correction in a batch script. Now I want to mosaic these 66 scenes which also basically works. Unfortunately at the left and right borders of the scenes there are very low values (here called 'no dada' values) which I cannot remove. There is a tool called 'S-1 Remove GRD Border Noise'. But applied it does not change anything.


texts: Import CSV to QGIS with z-values in column 4

texts:
I would like to import CSV into QGIS, but the z-values are listed in column 4. I didn't find out how I could specify that the values are not in column 3 but in column 4.

Edit: Added a screenshot


texts: API to build GoogleStreet-like photo maps

texts:
Do you know if Google offers a public API to display pictures in a map as they do with Panoranio?

Just to be clear, I don't want to use Google Street, or publish my pictures anywhere in the Internet, but just have blue dots as you can see in this map if you drag and drop the orange little guy.


texts: "AttributeError: getZ" error when trying to get the elevation of DXF linestrings

texts:
I'm using Python and OGR to extract the elevation of lines within a DXF file.

import ogr

driver = ogr.GetDriverByName('DXF')
datasource = driver.Open('test1.dxf', 0)

layers=datasource.ExecuteSQL( "SELECT DISTINCT Layer FROM entities" )
layer=datasource.GetLayerByIndex(0)

for i in range(0, layers.GetFeatureCount()):
        layerName = layers.GetFeature(i).GetFieldAsString(0)
        layer.SetAttributeFilter( "Layer='%s'" % layerName)
        print 'Layer=%s|Features=%s' % (layerName, layer.GetFeatureCount())

    for j in range(0, layer.GetFeatureCount()):
      myFeature = layer.GetNextFeature()
      myElevation = myFeature.GetGeometryRef().getZ(0)

But I get an "AttributeError: getZ" error when executing it.

What am I doing wrong?


texts: Python Toolbox Value Table values and list inputs

texts:
I am trying to create a Value Table I will use to create a field mapping object. First column will represent the target field, second column carries the source field. The first column is populated using the values property, so all fields from param0 are listed. I want the second column to offer a list of source fields from param2 to choose from.
I am able to see this list, but when I go to choose the matching field, I get:
The value is not a member of x. ERROR 000800: The value is not a member of CENTERLINEID | FROMLEFT | TOLEFT | FROMRIGHT | TORIGHT | FULLNAME |... etc.

Here is my code:
def getParameterInfo(self):
        """Define parameter definitions"""

    param0 = arcpy.Parameter('target_rcl', 'Target road centerline', 'Input', 'DEFeatureClass', 'Required')
    param0.filter.list = ["Polyline"]
    param1 = arcpy.Parameter('target_ap', 'Target address points', 'Input', 'DEFeatureClass', 'Required')
    param1.filter.list = ["Point"]
    param2 = arcpy.Parameter('rcl', 'Road centerline', 'Input', 'DEFeatureClass', 'Required')
    param2.filter.list = ["Polyline"]
    param3 = arcpy.Parameter('ap', 'Address points', 'Input', 'DEFeatureClass', 'Required')
    param3.filter.list = ["Point"]

    param4 = arcpy.Parameter('vt_rclfieldmap', 'Field mapping for RCL', 'Input', 'GPValueTable', 'Required')
    param4.columns = [['String', 'Target Field'], ['String', 'Source Field']]
    param4.filters[0].type = 'ValueList'
    param4.filters[0].list = ['x']
    param4.filters[1].type = 'ValueList'
    param4.filters[1].list = ['X']

    params = [param0, param1, param2, param3, param4]

    return params

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

    if params[0].value:
        flds = [f for f in arcpy.Describe(str(params[0].value)).fields if
                f.type not in ('OID', 'Geometry')]
        vallist = []
        for fld in flds:
            vallist.append([fld.name, ' '])

        params[4].values = vallist

    if params[2].value:
        srcflds = [sf.name for sf in arcpy.Describe(str(params[2].value)).fields if
                sf.type not in ('OID', 'Geometry')]
        params[4].filters[1].list = srcflds

    return

Here is an image that will hopefully demonstrate my goal... I feel like I'm close, just missing how to parse that a single value from the list of source fields is the desired input (edit) I believe the pyt takes my filter on the second column and applies it to the first, resulting in the error of 'value not found in list'. Not sure why it would disregard the index, but that seems to be the case.


texts: QGIS server don't show layers in Openlayers

texts:
I have a project with 3 layers , I see the layers names in the GetCapabilies and also as WMS layer in QGIS , But in Openlayers I get the tile empty.

I've used the following code :

var P4326 = new OpenLayers.Projection("EPSG:4326");
    var P900913 = new OpenLayers.Projection("EPSG:900913");
    var center = new OpenLayers.LonLat(34.87726678, 32.04670294).transform(P4326, P900913);
    var extent = new OpenLayers.Bounds(34.865,32.042,34.890,32.054).transform(P4326, P900913);
    var zoom = 15;
    var map, layer;

    function init() {

    var options = {
            controls: [

                    new OpenLayers.Control.Navigation(),
                    new OpenLayers.Control.ScaleLine(),
                    new OpenLayers.Control.PanZoomBar(),
                    new OpenLayers.Control.LayerSwitcher(),
                    new OpenLayers.Control.MousePosition()
                ],
        maxExtent: extent,
        restrictedExtent: extent,
        displayProjection: P4326
    };

    map = new OpenLayers.Map("map",options);
    var layerOptions = { zoomOffset: 11, resolutions: [76.4370282714844, 38.2185141357422, 19.1092570678711, 9.55462853393555, 4.77731426696777, 2.388657133483885, 1.194328566741943, 0.5971642833709715] };
    var base = new OpenLayers.Layer.OSM("Base", "", layerOptions);

     layer = new OpenLayers.Layer.WMS("WMS",
                "http://localhost:8081/qgis/qgis_mapserv.fcgi.exe",
                 {
                    layers: 'Water',
                    format: 'image/png',
                    map: 'D:/OSGeo4W/apps/qgis/bin/MyCity.qgs'
                 },
                 {
                    buffer: 1,
                    isBaseLayer: true,
                    graphicZIndex: 1
                 }
                );

    map.addLayers([base,layer]);

And only the OSM layer show , when I switch to the WMS , nothing shown.

Any Idea?!


texts: Azimuth greater than 360 in Veloview (Velodyne 32)?

texts: When I use Veloview, freezes at 1 of the frames and examine the data points, why are the azimuths greater than 360 degrees?


texts: Gdal won't read geotiff in parallel, but will read it in serial, and will read virtual raster in parallel. Why?

texts:
My algorithm takes a global dataset at 30m resolution, breaks it up into pieces, analyzes the pieces somehow, and writes the analyzed piece to disk. As all the pieces are independent, this is embarrassingly parallel.

I am using the joblib package to parallelize tasks. I have two global datasets: vrt.vrt (virtual raster) and tiff.tiff (geotiff). When I run a parallelized task that reads from vrt.vrt, I have no problems. However, when I run a parallelized task that reads from tiff.tiff, I get the following error message:

ERROR 1: LZWDecode:Wrong length of decoded string: data probably corrupted at scanline 368132
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: /path/landmasses.tif, band 1: IReadBlock failed at X offset 0, Y offset 368132
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 368132

The read command is the same for vrt.vrt and tiff.tiff:

tiff = gdal.Open('tiffpath')
start parallel
tiledata = tiff.ReadAsArray(
    xoff = int(xmin-npad),
    yoff = int(ymin-npad),
    xsize = int(ncol+2*npad),
    ysize = int(nrow+2*npad)
    )

However, when I run the tiff.tiff in serial, there are no errors and the output files look correct.

So does anyone know why gdal can't read from a .tiff using multiple threads, but can for a virtual raster?

EDIT:

In effort to progress, I built a virtual raster from the single geotiff file, and parallelization works fine. The question remains, however.


texts: QGIS - Layer missing polygons after pointsfile transformation

texts:
I'm using pointsfile transformation in QGIS for the georeferentiation of my layers. After the transformation the result some times has polygons missing like here:

Why is this happening? How do I prevent it?


texts: Bounding box styling of individual parts of multipart objects

texts:
In my QGIS project I have a multi-linestring vector layer where each multi-linestring feature consists of small not connected linestrings.

I would like to show a bounding box around each of the parts of each multi-linestring object like in the following example:

If I use the following styling

bounds($geometry)

I just get a bounding box around the entire multi-linestring object.

Is it possible to visualise the bounding box of each of the parts of the multi-linestring objects without creating a new single part dataset?


texts: How to get the direction raster from r.walk or r.cost

texts:
I am trying to use r.cost and r.walk and r.drain to calculate the shortest path from an archaeological site to the nearest basalt source. I am using all this algorithms with the processing toolbox (I have no grass experience...). I am able to run both r.cost and r.walk but when I try to use r.drain I have a python error. I thought I can use the rasters created to obtain the path in another program (Yes, like ArcMap). In order to do so I need the "movement direction raster" that you are suppose to get either from r.walk or r. cost. Nevertheless, when I use r.cost or r.walk from the processing toolbox and I can only visualize the cumulative costs rasters and not the direction ones. I also looked for an option to define a name and a path for the movement direction raster showing advanced parameters, but I have no luck either.
¿Does somebody know how to produce this movement direction raster? (if possible from the processing toolbox). I am using Dufour 64bits on win7 16mb of ram (standalone install and osGeo4w install)

Thanks in advance

Gerardo

PS. This is the Python error:
An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:\Users\IIAMaps/.qgis2/python/plugins\processing\gui\AlgorithmExecutionDialog.py", line 208, in accept
    self.setParamValues()
  File "C:\Users\IIAMaps/.qgis2/python/plugins\processing\gui\AlgorithmExecutionDialog.py", line 143, in setParamValues
    self.paramTable.valueItems[param.name]):
  File "C:\Users\IIAMaps/.qgis2/python/plugins\processing\gui\AlgorithmExecutionDialog.py", line 189, in setParamValue
    value.append(options[index])
IndexError: list index out of range

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.0.1-Dufour Dufour, d94c044

Python path: ['C:\Users\IIAMaps/.qgis2/python/plugins\processing', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python', 'C:\Users\IIAMaps/.qgis2/python', 'C:\Users\IIAMaps/.qgis2/python/plugins', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', 'C:\Program Files (x86)\ArcGIS\bin', 'C:\PROGRA~1\QGISDU~1\bin\python27.zip', 'C:\PROGRA~1\QGISDU~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISDU~1\bin', 'C:\PROGRA~1\QGISDU~1\apps\Python27', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', 'C:\Users\IIAMaps\.qgis2\python\plugins\mmqgis/forms', 'C:\PROGRA~1\QGISDU~1\apps\qgis\python\plugins\fTools\tools']


texts: OpenLayers: Switch layer off when switching another layer on

texts:
As the title says, I want to turn layers off when turning another layer on.
I wrote a function that works. However, I have to click the layer twice in the LayerSwitcher to turn the layer on that I want to see.

The first click turns the other layers off but only a second click turns the actual layer tho show on. Can somebody tell me where I am wrong in my code?

    function getVLayers(showLayer){
        for (var i=0; i < map.layers.length; i += 1){
            if (map.layers[i].visibility){
                if (map.layers[i].name == showLayer && map.layers[i].name == "OpenStreetMap"){

                } else {
                    map.layers[i].setVisibility(false);
                }
            }
        }
    };
    map.events.register("changelayer", this, function(e){
        getVLayers(e.layer.name);
    });


texts: Making Map with Color Ramp Using SQL Database Intergration?

texts:
This could be a trivial question, but I've been kind of stuck on it for awhile. I'm trying to construct multiple maps that reflect aggregated data showing overall risk by county in a SQL Database (using Microsoft SQL 2012). Each entry in my data set reflects an insurance contract with a risk score associated with a FIPS county code. I would like to use the AddJoin tool to join to an associated shape layer on FIPS county codes. The resulting map would show each county's risk with an associated color ramp depending on the risk level.

EDIT: I've taken RyanDalton's advice and created a seperate table that joins my data to a \ GIS database with shape entries on matching FIPS codes using MSSQL Server. I was able to connect to my specific database and make a new query layer useing ArcPy. The only problem I'm having is creating a color ramp based on the entries in my table (i.e. create a color ramp that shows differences between counties on values of Overall Risk Score).

The below image shows what I am trying to do with ArcPy, as well as the resulting map on the right (shows overall risk for an earthquake in the US).

Below is a rough sample of my code which doesn't output a map with a color ramp:

            import os, arcpy
            from arcpy import env

            # Declare Directories
            workingDir  = r'C:\Users\I55316\Documents\DQScore\ArcPy_Analysis'
            layerDir    = workingDir + r'\Layer_USNat'
            outputDir   = workingDir + r'\Output'
            dbDir       = workingDir + r'\DB'
            sdeFile     = 'DQ.sde'

            # Declare Database Server  Connection specifications
            databaseType        = r'SQL_SERVER'
            dataBaseInstance    = r'pmfldt43\sql2012'
            systemAuthType      = 'OPERATING_SYSTEM_AUTH'

            # Declare Query Layer Specifications
            PerilType               = 'EQ'
            layerName               = 'AllNationalCounty_EQ'
            query                   = 'SELECT * FROM dq.dbo.AllNational_EQ_GIS'
            oid_fields              = 'CNTY_FIPS'
            shape_type              = 'POLYGON'
            srid                    = '4326'
            spatial_reference       = arcpy.SpatialReference('WGS 1984 UTM Zone 12N')

            # Create Database connection via .sde file
            if os.path.isfile(dbDir+ '\\'+ sdeFile):
              os.remove(dbDir+ '\\'+ sdeFile)
            arcpy.CreateDatabaseConnection_management(dbDir,
                                                    sdeFile,
                                                    databaseType,
                                                    dataBaseInstance,
                                                    systemAuthType)

            arcpy.env.workspace = workingDir+'\\'+sdeFile

            arcpy.MakeQueryLayer_management(dbDir+'\\'+sdeFile,
                                          layerName,
                                          query,
                                          oid_fields,
                                          shape_type,
                                          srid,
                                          spatial_reference)

            # Declare map document
            mxd = arcpy.mapping.MapDocument(workingDir+r'\NationalDQScoring_EQ.mxd')
            mxd.author = "Krishna Soni"

            # Add Layer
            #lyrCounties    = arcpy.mapping.Layer(layerDir + r'\UnitedStates_Counties_TS2.0.0.lyr')
            lyrFaultLines  = arcpy.mapping.Layer(layerDir + r'\UnitedStates_FaultLines_TS2.0.0.lyr')

            # Define data frame

            df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]

            # Add Layers
            arcpy.mapping.AddLayer(df,lyrFaultLines)

            # Add graduated Colors
            for lyr in arcpy.mapping.ListLayers(mxd):
              if lyr.symbologyType == 'GRADUATED_COLORS':
                  lyr.symbology.valueField = 'Overall_Risk_Avg'
                  lyr.symbology.numClassed = 10


texts: How to display faster shapes covering a country map?

texts: I need to display ~9000 polygons that represent country regions on a google map in a fast way. The polygons are kml formats (but I can get other formats too). At the moment everytime the map is displayed, the relative shapes are fetched out from a database. However there is a sort of delay which is not optimal. Is there anyway to avoid this database fetching or some optimisation maybe fully frontend side?


texts: Using QGIS and Python to develop sound attenuation model?

texts:
I am in the process of developing a model that allows me to see sound attenuation over distance from a point source. I am using this sound attenuation equation:

Lp = Lw - 20logr - 11 - Aabs - Ae

where Lp = sound pressure (dB); Lw = sound pressure at point source (db); r = radius from source (m); Aabs = atmospheric absorption; Ae = excess attenuation (vegetation, barriers...)

I have a DSM over flat ground that I would like to apply this model to, it has open grassland, vegetation and houses. I will use conditionals to insert the appropriate variable for the surface type.

Is there a simple way to apply the equation to each pixel in a raster from the point source, iterating out radially (not spirally) instead of in a loop? As I am using the DSM, the newly generated value is directly dependent on the previous db(A) value in the line of sight from source. I should create a new value for each pixel until I reach background noise levels - this new data can be used to create and overlay a noise map on to a basemap.

If the approach seems way off, please be critical -- up until now I assumed that I should be looking for a way to apply this equation to a NumPy array; perhaps this is also an inefficient methodology? (Python newbie!)


texts: ArcGIS 10.2.2 Chart Graph data precision shown

texts:
I created a pie chart that had several pieces of data; pipe diameter to length of pipe per diameter. I summarized data from length of pipe and diameter from a database and input into a summary data table. I then tried to graph the table and the data shows up with 3 decimal points for the summarized lengths. I modified the editor options to indicate no decimal points and changed the table setting to show no decimal points.

Is there a simple way to change the precision of the data in the graph?


texts:
How would I go about labeling every other feature in QGIS (1.8)?

For example I have a set of contours every 10 meters from 0m to 100m. How would I label only every 20 meter contour? e.g. label contours 0, 20, 40, 60, 80 and 100m while showing all the contours.


texts: i am working on grass gis. how to use grass gis for liss4 of IRS-P6.I have to do atmospheric correction for hilly terrain data.i would like to use grass's i.atcorr. may i know the dependencies and how to customize it for my problem.will liss4 data directly used in that module?or only it is limited to specific sensors. i would also like to inform that i have compiled OTB4.0 but the code assumes that surface as Lambert which not suitable to my study. hence i want to use grass gis for it. will grass also has same assumption of lambertian surface? please suggest me.


texts: Updating PostGIS layer from another vector layer in QGIS?

texts:
I have a vector layer that contains way points that I load using "load vector layer" from a gpx file.   I am attempting to add these new points to an existing postgis layer by simply copying and pasting between the layers.

When I first tried this I got an error "Column has Z dimension but geometry does not".  Sure enough the PostGIS table was 2D so I converted the geometry to pointz and retried but I still get the same error?

database manager now shows the the_geom of the postgis table as

6  the_geom geometry (PointZ,2193)

How can I inspect the geometry of a layer from within QGIS? Any other suggestions on how to figure out what is really wrong.


texts:
I need to project data which can be nearby the poles, and I want to use a custom transverse mercator projection for a number of reasons.

Using GeographicLib (either TransverseMercator or TransverseMercatorExact), I see I can only set a projection origin for longitude, but not for latitude. However, what kind of projection accuracy can I expect at the poles?

/**
 * Forward projection, from geographic to transverse Mercator.
 *
 * @param[in] lon0 central meridian of the projection (degrees).
 * @param[in] lat latitude of point (degrees).
 * @param[in] lon longitude of point (degrees).
 * @param[out] x easting of point (meters).
 * @param[out] y northing of point (meters).
 * @param[out] gamma meridian convergence at point (degrees).
 * @param[out] k scale of projection at point.
 *
 * No false easting or northing is added. \e lat should be in the range
 * [&minus;90&deg;, 90&deg;].
 **********************************************************************/
void Forward(real lon0, real lat, real lon,
             real& x, real& y, real& gamma, real& k) const;


texts: OpenLayers: Copy shape layer in a vector layer

texts:
I have a web-gis application based on OpenLayers and GeoServer.
I need to modify a shape layer. I'm able to show the shape using the OpenLayers.Layer.WMS method. In fact I published the shape in GeoServer and called it in the OpenLayers map file.

My idea is to copy the shape layer in a vector layer. In this way I can modify the vector layer before to save it.

How can I copy the shape layer in a vector layer?


texts: Postgres/PostGIS skip intersect on features already intersected

texts:
I have a collection of polygons(building outlines called "Build_poly") that I am intersecting with another table of flood defenses (also polygons called Flood_poly).
When ever one of the building polygons intersect with a flood defense polygon I am updating the "defense" field of building polygon from 0 to 1.

What I have been trying to do to speed up my query is put in a condition where: if a building_poly has successfully intersected with a flood polygon for it to skip trying to intersect with any other flood defense polygons. This is what I've got so far, but I haven't managed to get any form of skip to work as this is my first trying to use postgresql.

UPDATE Build_poly
SET Defenses = 1
FROM Flood_poly
WHERE ST_Intersects(Build_poly.geom,Flood_poly.geom)"


texts: Identifying and Placing Points on DEM Features using ArcGIS Desktop?

texts:
I am using UAV acquired imagery to identify trees at a peat land in the UK. I have a DEM and an Orthomosaic which i have created which can bee seen here  (note: the Orthomosaic is at a much higher resolution when viewed in ArcMap than online).

Manually i can identify trees by looking at the high points on the DEM and by visually spotting them on the Orthomosaic and placing a point on them.

How would i start going about getting ArcMap to do this automatically for me?


texts: Calculate BBOX by latitude and longitude

texts:
How openlayers3 generate BBOX by latitude and longitude.

I want to generate BBOX for GetFeatureInfo request by latitude and longitude in C#.


texts: Google Earth Engine - Error downloading cloudfree Sentinel image

texts:
I got the following error message after trying to import cloud-free Sentinel-2 imagery for a polygon area covering the entire Tiwi Islands, Australia.

  true-colour image: Layer error: internal error

// Filter to only include images intersecting Tiwi.

var polygon = ee.Geometry.Polygon({
  coords: [[
[129.90415998438357,-12.027354108834277], [131.63999982813357,-12.027354108834277],
[131.63999982813357,-11.101721188177093], [129.90415998438357,-11.101721188177093],
[129.90415998438357,-12.027354108834277]]],
  geodesic: false
});

//  Define the image collection we are working with by writing this command
    var image = ee.Image(sent2

    // filter to get only images in the date range we are interested in
    .filterDate("2018-09-01", "2018-11-30")

    // sort the collection by a metadata property, in our case cloud cover is a very useful one
    .sort("CLOUD_COVERAGE_ASSESSMENT")

    // select the first image out of this collection - i.e. the most cloud free image in the date range
    .first());

    // print the image to the console.
    print("A Sentinel-2 scene:", image);

    // Define visualization parameters in a JavaScript dictionary for true colour rendering. Bands 4,3 and 2 needed for RGB.
    var trueColour = {
        bands: ["B4", "B3", "B2"],
        min: 0,
        max: 3000
        };

  // Add the image to the map, using the visualization parameters.
  Map.addLayer(image, trueColour, "true-colour image");


texts:
Is it possible to create a buffer polygon using an attribute field which is in different units than the feature class coordinate system?

Example: I have a field that stores distances in miles (double data type) and the feature class is in state plane (units are feet).

Is there a way to do this?

I'm using Esri ArcGIS 10.2.


texts: Extracting climate data for given coordinates from interface?

texts:
I have my own dataset of latitudes and longitudes. I need to extract annual average climatic variables (Rainfall, temperature,sunshine etc) for this specific data set. I tried WorldClim but it requires me to download other software to process it. I tried to use these databases from Is there a simple interface from which I can extract average summer temperature and humidity for a list of latitude and longitudes? but it appears that I need to input each value of long-lat to get climate data here at Weather Underground.

Is there an easier database, where I simply import my lat-longs from excel/csv and get the required climate data?


texts: Can I overlay layers in Dealul Piscului 1970 and WGS84 in QGIS?

texts:
I have a 10 Gb ECW file (mosaic of a whole county) that is in Oblique Stereographic projection with Datum: Dealul Piscului 1970. The problem was that no matter the Default CRS for new project was, when I wanted to add this raster the coordinates were all wrong, although the input projection for the raster was correct. This problem was solved by deleting the ECW's ".xml" file that was accompanying the ECW file. Once deleted, the problem was solved and the QGIS was showing the right coordinates.

Once the above problem was solved. I wanted to add a vector layer of 30 mb that represents the whole waterway of the Danube Delta plus inland towns. This layer is in WGS84. Now when I try to import this layer does not overlay with the ECW file. The ecw file have been defined with the right projection also the vector layer have been defined with the right projection. In settings/options/CRS/ the "on the fly" is checked and the project projection is WGS84 and also tried the vice-verse with Oblique Sterographic Dealul Piscului. With all the permutations available I still can not overlay the layers.

Do you have any other suggestions?

Thank you


texts:
I am trying to calculate the amount of water in this reservoir. The green dots are water depths in meters (possitive value). The water depth at the edges of the reservoir is 0.

How can this be done using QGIS?


texts: Getting libpng error for QGIS on Ubuntu 19.04?

texts:
So I have already attempted the solutions that were posted here, and I still get the same issue even after reinstalling QGIS: libpng problems for QGIS

I keep getting the following error message:

libpng warning: Application was compiled with png.h from libpng-1.6.36
libpng warning: Application  is  running with png.c from libpng-1.2.56
libpng error: Incompatible libpng version in application and library

And then when I try to open anything with QGIS, it crashes and then I get this message in the terminal:

libpng warning: Application jmp_buf size changed
QGIS died on signal 11Could not attach to process.  If your uid matches the uid of the target
process, check the setting of /proc/sys/kernel/yama/ptrace_scope, or try
again as the root user.  For more details, see /etc/sysctl.d/10-ptrace.conf
ptrace: Opération non permise.
No thread selected
No stack.
gdb returned 256
Abandon (core dumped)

I have tried to rebuild GDAL, which I have installed from source, with updated libpng headers, but that didn't fix my issues. I also installed QGIS from the following apt repository:

Atteint :8 https://qgis.org/ubuntu disco InRelease

And with my installed QGIS, I tried the following command, and I got nothing back:

ldd `which qgis` | grep png

I am at a loss as to what to do, short of completely reininstalling Ubuntu.
 Has anyone else encountered this issue and found a way to fix it?

EDIT: This is not a duplicate question because there are additional error messages that came with this particular issue, and the fixes that addressed the previous so-called "duplicate" question, as I stated before, did not work. If they worked, and if there weren't also new error messages that make this situation different, then I wouldn't be asking it.


texts: Python and postgis, trying to figure out data types

texts:
I have a simple sql I am running from python that gets all these items from a table (mostly the mpoly I am curious about).  Which is a geometry.

cur = self.db.cursor()
cur.execute('select St_AsText(ST_Centroid(mpoly)) from swsite_cesiumentity')
allEntities = cur.fetchall()

Then I can get the info out:

    print ' ***'
    print allEntities[0]
    print '----'

and it will give me:

Row(st_astext='POINT(-87.266150008217 30.0857181963441)')

That is well and good, but I try to use this data to then get all the points that are so far away from the above point and I do not know how to use the types I have in the sql query..  What I thought should work was:

for currentEntity in allEntities:
    print currentEntity
    cur2  = self.db.cursor()
    cur2.execute ('SELECT * FROM swsite_cesiumentity WHERE ST_Distance_Sphere(mpoly, ST_GeomFromText(' + currentEntity.st_astext + ')) <= radius_mi * 1609.34)) from swsite_cesiumentity')
    list = cur2.fetchall()
    print list

It doesn't like what I am feeding the second sql statement:

pg.ProgrammingError: ERROR:  syntax error at or near "30.0857181963441"
LINE 1: ...ere(mpoly, ST_GeomFromText(POINT(-87.266150008217 30.0857181...
                                                             ^

C02RH2U9G8WM:ingest sth$


texts:
I am looking to create a PDF using Arc10 that contains locations for our site. At each of these locations you would be able to click the point and you can then see the photos for the area. Is this possible?

I am thinking something like Google Earth's functionality for shared photos. I feel like this is feasible.

Thanks!

Tony


texts: Remove vertices coincident with a point layer

texts:
Is there any way to delete vertices from a feature which are coincident with a point layer? I mean in an automatic way.

I have many lines forming a "T shape", so I would like to remove the middle vertex of the horizontal line so the only vertex there is the first (or last) one from the vertical line. Every line has a different lenght, so it's not possible to smooth or simplify the feature.

I have to do it with QGIS (but if there is a way with another software, i would also appreciate the answer).


texts: Buffer in meters instead of degrees in south america

texts:
I need to create some buffer in meters, but it not possible how I'm doing. I really don't know what is happening.

My code:

library(rgdal)

minas_ti<- readOGR("minas_TI.shp")

minas_ti<-spTransform(minas_ti, CRSobj = "+proj=aea +lat_1=-2 +lat_2=-22 +lat_0=-12 +lon_0=-54 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")
proj4string(minas_ti)

buff_min_ti10km<- gBuffer(minas_ti, width=10000, byid = T)

When I use the code above the buffer is generated in degrees.


texts: Styling polygon layer with symbols that indicates value of attribute of polygon?

texts:
I want the same thing, but the number of dollars shown should indicate the value of the attribute of the polygon. For example one polygon may have 2 dollar icons indicating 200 (value of an attribute, say income) and another polygon shows 5 icons indicating a value of 500.

The symbols should be automatically generated indicating the rounded value of the attribute, and they should be in a single line.
I tried using the centroid fill with single marker, but could not understand which setting I need to play with to achieve this effect.


texts: Loading Json properties in a popup using leafletjs

texts:
I have a geojson file (tide_gauges.geojson)

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "NUMERO": "10569", "MAREGRAFO": "COCALZINHO", "FICHA": null }, "geometry": { "type": "Point", "coordinates": [ -4, -1 ] } },
{ "type": "Feature", "properties": { "NUMERO": "50240", "MAREGRAFO": "PONTA DE PARANAPUA", "FICHA": null }, "geometry": { "type": "Point", "coordinates": [ -4, -2 ] } }...

I retrieved this Json file in a variable:

var estacoes = new L.geoJson();
estacoes.addTo(map);

$.ajax({
dataType: "json",
url: "geojson/tide_gauges.geojson",
success: function(data) {
    $(data.features).each(function(key, data) {
        estacoes.addData(data);
    });
}
}).error(function() {});

I am in need to set all the feature properties (NUMERO,MAREGRAFO,FICHA,LATITUDE and LONGITUDE) in a popup that appears after a click over point feature on the map.

I tried the code below, but it is not enough:

var popupContent ="";
for (var k in feature.properties) {
    var v = String(feature.properties[k]);
    popupContent += k + '->' + v ;

Could someone give some tips on that, please?


texts: How to export/save a layer composition to GeoTiff?

texts: I want to export/save a QGIS composition of raster/vector layers to a GEOTIFF or ECW raster with high resolution (20 000 * 20 000 Pixel).


texts: Selecting features in QGIS via NetworkAPI HTTP POST request

texts:
Searching for a way to remotely operate on a running QGIS instance I recently stumbled upon a GSOC 2017 project NetworkAPI which provides a REST interface for QGIS (https://gitlab.com/qgisapi/networkapi). Installation and creating a new project via http request works fine, directly from Firefox adress bar as well as with Firefox REST Easy addon:

http://localhost:8090/qgis/newProject

What I want to do is remotely select features of a certain layer but unfortunatly NetworkAPI does not implement mapLayer.selectByExpression() directly. Thus I tried to use

http://localhost:8090/qgis/processing/runalg?alg=qgis:selectbyexpression

with REST Easy at first with no parameters, which returns an error

[
  "Error: Wrong number of parameters",
  "ALGORITHM: Select by expression",
  "\tLAYERNAME <ParameterVector>",
  "\tEXPRESSION <ParameterString>",
  "\tMETHOD <ParameterSelection>",
  "",
  "",
  "METHOD(Modify current selection by)",
  "\t0 - creating new selection",
  "\t1 - adding to current selection",
  "\t2 - removing from current selection",
  "\t3 - selecting within current selection",
  ""
]

as expected which basically means to me that this approach is at least promising.

Following the documentation from line 101 of https://gitlab.com/qgisapi/networkapi/blob/master/functions_processing.py

"...The arguments to the processing algorithm should be passed in JSON format in the POST body of the request, and the request MUST specify 'Content-Type: application/json' in its header. The body should either be a simple JSON array (passed to the algorithm as positional arguments in that order) or as a JSON object containing all desired arguments as 'name: value' pairs...."

I tried to set up REST Easy as shown in the image below

but this is just causing a NetworkAPI connection timed out afer 3000ms.

NetworkAPI Log Messages panel:

2018-09-10T16:02:34 0   Listening on port 8090
2018-09-10T16:02:34 0   Processing connection #6 (127.0.0.1)
2018-09-10T16:02:34 0   Connection opened...
2018-09-10T16:02:34 0   readyRead() was signalled before full HTTP request line was available for reading (0 bytes in buffer)
2018-09-10T16:02:34 0   Still busy with #6, putting incoming connection on hold...
2018-09-10T16:02:34 0   Connection opened...
2018-09-10T16:02:34 0   Processing request: received 3 of 423
2018-09-10T16:02:37 1   Connection timed out after 3000ms
2018-09-10T16:02:37 0   Disconnecting #6 (127.0.0.1)
2018-09-10T16:02:37 0   Listening on port 8090

Do I do the POST request the right way and are my parameters correct (tried with and without quotes, no difference)?


texts: Convert a SHP file set to a categorised ASC raster with Qgis?

texts:
I have a shapefile with regolith (soil) data that I would like to use as a categorised raster (say 007_1ASC) in Maxent.

The regolith data is in the following files:

007_1.cpg
007_1.dbf
007_1.prj
007_1.sbn
007_1.sbx
007_1.shp
007_1.shx

Using Qgis, how can I get this into an ASC file with categorised data?


texts:
I would like to show the edges of the result route in a table, like this:

http://routing.geosdi.org/

I did the pgrouting part, just have to make the table and connect it.

How can I achieve this?


texts: Automatically fill table with Geographic coordinates if shapefile is in UTM?

texts:
I have a shapefile of points in Qgis 2.18. I created the points by direct click on a map, using the digitizing tool. Now I need to fill the 2 columns "Lat" and "Long" in the attribute table with Geographic coordinates (DD MM.mmm).

The shapefile is in Utm. I used the instructions found here: How to automatically fill a lat/long in attribute table.

But when I try to use the Field calculator, the columns are filled with Lat and Long values referred to UTM. Even when I set the Project properties to Geographic, the Field calculator calculate always UTM coordinate to fill the columns.


texts: Raster not showing on Orthographic projection when zooming out

texts:
I am trying to put together Topographic .Tif rasters (from GMTED data) on an ortographic globe on Qgis 3.4.1, for exemple these 8 layers:

However when I change the projection and one of the various layer appears toward the extremities of the globe, it suddenly dissapears: (the above screenshot used a projection centered on 35°W, the one below on 60°W)

However If I zoom/increase the scale on the rasters that dissapeared, they eventually reappear

It doesn't matter whether I change the scale or zoom, as long as I unzoom at a certain level it disappears.

I encounter the same problem when I make a print layout and add a part of the map on it, if it is unzoomed then the raster disappear.

I tried to check in "scale dependent visibility" in the render part of the layer settings, but it didn't do anything. I also don't think it is a problem with parts of the raster being "behind" the globe and thus not being rendered, since I encounter this problem even when all corners of the raster are on the visible half of the globe, although this glitch happens nearly always on rasters that are near peripherical to the globe.

Thanks for your help.


texts: Import of OSM data into PostgreSQL using osm2pgsql: which formats are supported?

texts:
I followed a tutorial to install a PostgreSQL server on my Mac OS X Mavericks, and to import some OSM data with the osm2pgsql converter (https://www.mapbox.com/tilemill/docs/guides/osm-bright-mac-quickstart/).

I used OSM data in the .osm.pbf format for the import. That works fine.

Is it possible to use the formats .osm and .osm.bz2 as well? And maybe other formats?
Do I then need to add certain options to the osm2pgsql command, depending on the chosen format?

I intend to import only small amounts of data into the database, mainly for visualizing purposes in TileMill.


texts:
I am at the mercy of a govt authority who have provided utility information in a cad format with a 0,0 origin and a nonsense grid. About a dozen files, the grid and origin are entirely arbitrary.

To do an affine on a vector layer in QGIS, I propose to the use the Geoprocessing Plugin - Vector Affine Transformation.

But my math skills betray me....

The matching points for the ground control are

pt1 raw 6.42775 3.61140
pt1 gcp 567742 6450622

pt2 raw 2.7524 4.1539
pt2 gcp 562285 6450858

pt3 raw 6.0293 6.4965
pt3 gcp 566709 6454802

(edited x 2 , there are now 3 points here) A little rough on the correct locations.

The ground control is a UTM format.

What would be the form of the equation, so as I can set it up in a spreadsheet and plod away?

The input for the the plugin is Scale X Scale Y Rotation X Y (not used I expect) Translation X Translation Y


texts: IndentationError: expected an indented block

texts:
I am getting the IndentationError in arcpy. I understand this is covered elsewhere but none of the tips appear to be working, I cannot figure out the error.

Here is my code;

def getCoords(self):

    with arcpy.da.SearchCursor(self.inFC, ['SHAPE@XY']) as rows:
        self.coordList = [r[0] for r in rows]

I am getting the error at the line that starts "with". I am using Notepad++ and have changed tabs to spaces (4) in preferences. I have deleted and re-written my code a number of times to ensure no erroneous spaces, however I am still getting this error.

self.coordList is created when the class is initialised;

def __init__(self, inFC, buffDist, sector): #Must have a set of points and a value to buffer these by. Also, sector.

    self.inFC = inFC
    self.buffDist = buffDist
    self.sector = sector
    self.coordList = []

What mistake am I making?


texts: Weighted Average of several raster "zones"

texts:
Suppose I have two rasters, A and B. I want to determine the value of a certain "zone" in raster A by taking the weighted average of several other zones from raster B, that partially overlay onto that zone of raster A.

I think its much better if I show you a picture to try to explain:

The "certain 'zone' in raster A" is shown as the teal polygon in the border, suppose it is rasterized and is made up of 1000 cells.
The red grids are the "several other zones from raster B." In this case, there are four denoted by g1, g2, g3, g4 (these are some values).
The weight will be determined by the number of pixels lying within the respective zones in raster B.
The value of that certain 'zone' in raster A will be determined by this:

value of zone in A = (50/1000)g1 + (300/1000)g2 + (150/1000)g3 + (500/1000)g4

I've thought of employing the Combine tool in ArcMap and then bringing it into excel.. but it is a fairly manual process. I have a lot of data to get through (like 50+ raster B's that I must aggregate through)..

Is there an easier way to do this? Is there a specific function you think may help?
I have little experience in writing a python script but if you point me to some packages+functions in arcpy or numpy or whatnot, I'd greatly appreciate it.

====

edit: Similar, previous questions have suggested things like "dasymmetric mapping" like this thread: Using National Agriculture Statistics Service (NASS) county level crop data to estimate crop acreage in hydrologic regions?.. which suggests the use of the Intersect tool, but I'd generally prefer to rely on using weights from the counts of defined cell sizes of a raster rather than have arcgis calculate the area proportions from vector data.

I also found a helpful tip on this question, Obtain normalized raster average per polygon?, suggesting Zonal Statistics. Is the "mean" of Zonal Statistics the same as what I'm trying to do? I'll try it and report back...


texts: Distance matrix (nearest neighbour) with polygons

texts:
I would like to calculate the closest distance between polygons in the same layer. The "distance matrix" tool only accept point layers.

Is there a (easy) way to calculate the distance between edges and edges ?

I work with qgis 3.2


texts:
I'm currently trying to cutline crop a tif with a shapefile and no matter what I do I have pixel either overhanging the cutline or being cropped too early.

So here is an example picture of what is happening when I use gdalwarp:

When I use arcpy to crop I get this result:

The code I am currently using with gdal warp is:

 gdalwarp.exe -ot "float32" -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84
 +no_defs " -t_srs "+proj=utm
 +zone=48 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "-of "GTiFF" -cutline
 mask.shp -crop_to_cutline infile.asc outfile.tif -wm 5000 -tap -tr 100 100 -overwrite

I am quite keen to use gdal if possible to do this cropping as it is considerably faster

I have been asked to provide gdalinfo info, see below.
This is from the gdal code:

This is from the arcpy code:

SOLUTION was found in comments.

Exactly what I wanted:

My new problem to ponder next is:
[]


texts:
I want to get attributes of feature using getFeatureInfo:

if(sss==null){
   if (sld==null){
     sld= new OpenLayers.Control.WMSGetFeatureInfo({
    url: 'http://localhost:8080/geoserver/wms',
    title: 'Identify',
    layers: myLayer,
    queryVisible: true
     });
   app.mapPanel.map.addControl(sld);
   sss=app.mapPanel.map.getControlsBy('title', 'Identify');
}
}
for(var i=0;i<sss.length;i++){
alert("BROMACIA");
 sss[i].activate();
}
sss[0].events.register("getfeatureinfo", app.mapPanel.map, function showInfo(evt) {
alert(evt.text)});

And nothing heppens.
When i try to do in firebug:

>>> sss[0].events.register("getfeatureinfo", app.mapPanel.map, function showInfo(evt) { alert(evt.text)});

firebug says undifined.

Whats wrong?


texts:
I try to show a small resolution raster file in the QGIS Print Manager tool.
I would like to show a small label with each cell value (position - circa middle of each cell).

How I can do it? Is there any plugin which can help me to accomplish this task?

My idea is to use gdal2xyz, load such file using CSV plugin and then save it as a vector file.  Unfortunately it is much work, because I have many raster files.

Another idea is to use gdal2xyz.py, then ogr2ogr to save it as vector file.

Is there any way to batch this process? Is there any open-source program where I can do it much more easily?


texts: QGIS 1.8 Vector layer from ODBC connection

texts:
I have successfully set up an ODBC link from a FileMaker DB to QGIS, and using Layer/Add Vector Layer I can get a layer with a fully populated attribute table listed in the Layers Panel.

I now want to plot datapoints from the new layer but cannot work out how to tell QGIS that the new layer has georeferenced information and which fields contain the LAT/LONG data.  I guess that I am expecting to find a window similar to that which appears when you create a layer from a Tab Delimited file, where you define which fields contain the LAT/LONG data.

Can anyone help me out?

Thanks


texts: How to get feature information from a WFS at a certain coordinate (in Python)?

texts:
I have a WebFeatureService at a remote url. It contains different 'layers'/ types.
How do I get information at a certain point (lat, lng) similar to what QGIS returns when you click the map using the 'Identify Features' option (see screenshot)?

The getfeature method of owslib only lets me input bounding boxes (bbox) and not a coordinate. Is there a simple way to get this information?


texts: How to update a popup content from its marker?

texts:
I created a popup with Leaflet doing so:

marker.bindPopup(content).openPopup();

How can I update the content value afterwards?

I suppose to do it from the marker, something like this:

marker.updatePopup(newContent);


texts:
I am looking for help. I have done a little bit of geo-registration of images with QGIS. But now, I have been handed some LIDAR data in PCD format (basically x,y,z), But it is not geo-referenced.

I have looked at: Georeferencing LAS files in geographic coordinates (lat/lon) or in cartesian (projected) ones?, but I do not have LAS data, and that article talks nothing about what would work for a GUI to do GeoReferencing, as QGIS provides.

I know the exact area and have obtained a GeoTiff from Google Earth. The entire region is only about 30 acres. So I am assuming that a simple linear transform of the values will be sufficient. But I am was wondering if any GUI tools could help with the conversion (at least getting the ground control points to line up and the linear formula). I have used QGIS to do this for unknown raster data, but this is the first time I have been asked to do this for point cloud.

My simple research has found a ideas, which likely won't work, but maybe someone can tell me what does work:

I have found an interesting free tool called cloud compare which if I understand the demos, I might be able to use (it looks like it needs two point clouds, but I could convert the GeoTiff with FME to a point cloud). I do not have ESRI available. Is Cloud Compare able to do this? How hard would it be?

If we want to put out $400 perhaps the Blue Marble folks can do this with their LIDAR plug in: http://forum.globalmapperforum.com/discussion/7444/georeference-static-lidar-data , does anyone recommend that?

I have also see this in the Matlab Vision Toolkit: http://www.mathworks.com/help/vision/ref/pcregrigid.html#bux42c5


texts:
I'm trying to adapt the pgr_fromAtoB function to pgr_fromAtoMult replacing pgr_dijkstra function to pgr_drivingdistance.Here is the code:

CREATE OR REPLACE FUNCTION pgr_fromAtoMult(IN tbl character varying, IN x double precision, IN y double precision, OUT seq integer, OUT gid integer, OUT name text, OUT heading double precision, OUT cost double precision, OUT geom geometry)

  RETURNS SETOF record AS
$BODY$
DECLARE
    sql     text;
    rec     record;
    source  integer;
    point   integer;
    dist    double precision;

BEGIN
        -- Find nearest node
        EXECUTE 'SELECT id::integer FROM layer_vertices_pgr
                        ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
                        || x || ' ' || y || ')'',4326) LIMIT 1' INTO rec;
        source := rec.id;

        -- Changing pgr_dijkstra to pgr_drivingdistance query
        seq := 0;
        sql := 'SELECT id, geom_way, osm_name, cost, source, target,
                                ST_Reverse(geom_way) AS flip_geom FROM

' ||
                        'pgr_drivingdistance(''SELECT id as id, source::int, dist::double precision, '
                                        || 'time_4::float AS cost FROM '
                                        || quote_ident(tbl) || ''', '
                                        || source || ', ' || dist
                                        || ' , false, false), '
                                || quote_ident(tbl) || ' WHERE id2 = id ORDER BY seq';

        -- Remember start point
    point := source;

            -- Flip geometry (if required)
            IF ( point != rec.source ) THEN
                    rec.geom_way := rec.flip_geom;
                    point := rec.source;
            ELSE
                    point := rec.target;
            END IF;

            -- Calculate heading (simplified)
            EXECUTE 'SELECT degrees( ST_Azimuth(
                            ST_StartPoint(''' || rec.geom_way::text || '''),
                            ST_EndPoint(''' || rec.geom_way::text || ''') ) )'
                    INTO heading;

            -- Return record
            seq     := seq + 1;
            gid     := rec.id;
            name    := rec.osm_name;
            cost    := rec.cost;
            geom    := rec.geom_way;
            RETURN NEXT;

    RETURN;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100
  ROWS 1000;
ALTER FUNCTION pgr_fromAtoMult(character varying, double precision, double precision)
  OWNER TO postgres;

And here is the Error:

ERROR:  el registro «rec» no tiene un campo «source»
CONTEXT:  sentencia SQL: «SELECT ( point != rec.source )»
función PL/pgSQL pgr_fromAtoMult(character varying,double precision,double precision) en la línea 32 en IF
********** Error **********

ERROR: el registro «rec» no tiene un campo «source»
SQL state: 42703
Context: sentencia SQL: «SELECT ( point != rec.source )»


texts: QGIS modelbuilder, having a choicebox/value map of strings as string input

texts: I was wondering how you can get a choicebox/value map in the string input? Right now we have to type in value in but we would like to have it that has a choicebox/listbox to save time with typing.


texts: SQL query fails in arcpy, but works in ArcGIS Pro 'select by attributes'

texts:
I'm building an SQL query in Python that is used to make a feature layer, but the feature layer returns no features, when it should return 1169 features.  The query below works on 95% of the data (1955-2017), but fails to find the 2018-2019 data.  When I copy the exact string printed in the python console into the SQL 'select by attributes' window in ArcGIS Pro, it works as expected (see screenshot).  Code sample is:

for year in years:
    for event in types:
        where_clause = f"Year = '{year}' and Type = '{event}'"
        print(where_clause)
        arcpy.MakeFeatureLayer_management(spc, "temp_lyr", where_clause)
        print(f"Number of {event} events in {year}: {arcpy.management.GetCount('temp_lyr')}")

for example, returns:

Year = '2018' and Type = 'Tornado'
Number of Tornado events in 2018: 0

The original data (1955-2017) was imported into a feature class from a shapefile. The 2018-2019 data was added with an insert cursor from a CSV that had lat/lon information and the same fields.  That seems to be the only difference in rows that work and don't work in arcpy.  They are text fields and I've already cleaned up all whitespace.  Any ideas what could be going wrong?

The data is in a file geodatabase.  I'm using ArcGIS Pro 2.4.2 with Python 3.6.6.

p.s. I also already tried wrapping the query in a docstring and adding double quotes around the field name.

where_clause = f""""Year" = '{year}' and "Type" = '{event}'"""

which returns something like:

"Year" = '2018' and "Type" = 'Tornado'

Here's a sample of the table.  In my arcpy script the first two rows (2017, Tornado) will be included in the layer in the correct iteration of the loop.  The next two rows (2018, Tornado) will not be included.  All fields are text except for Lon/Lat.


texts: Algorithm to find the polygon ID where a point belongs to

texts:
I am wondering what is the simplest way to perform a look up of polygon Id based on the point location.

Input: a point with coordinates, a table of polygons

Output: the polygon Id where a point is located.

I am thinking to preload the polygon table into PostGIS (The size is about 220,000 entries). What is the best way to process all these data?

Thanks.


texts: What is happening different to my TIFF between export data and arcpy

texts:
I am working in FME with TIFFs downloaded from sentinel-1 (SAR). These TIFFs display correctly in ArcMap but in FME is loses its spatial reference and a bound box calculation gives the number of pixels in the tiff from the upper left corner at 0,0.

I concluded that this is an issue with the flavour of TIFF, so I exported the tiff in ArcGIS to a new file. FME could read this fine. A bounding box output would appear in the right place.

As I'm trying to automate a workflow, a manual export in arcmap won't be good enough. I've tried using the Python Caller in FME with

def processFeature(feature):
    origraster = feature.getAttribute('path_windows')
    print "!!!!!!!!!!!!!!!!!!!!!!!!!!! RECREATING RASTER!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
    print origraster
    newraster = origraster.replace(".tiff", "_arcgis.img")
    print newraster
    arcpy.CopyRaster_management(arcpy.sa.Raster(origraster),newraster)
    feature.setAttribute("newraster", newraster)

    pass

This creates a new raster (I've tried .tif, this is another test of using ERDAS Imagine .img) but one with has the same problems as the original.

I've tried RasterToOtherFormat_conversion but have the same problem.

Is there a difference between the process in exporting from ArcMap and copying in arcpy? And if so how can I remove this?

This is the info given in FME data inspector of the Python produced .img

And this is what it sees if the .img is produced by manually exporting in ArcMap


texts: Setting ArcMap attribute dialog split orientation to horizontal by default?

texts:
When you open ArcMap 10, and open the attribute window (editor toolbar), the attribute dialog is split vertically.

You can press the toggle button to split the dialog horizontally.

I always use it split horizontally, so how can I set it to that by default?


texts: Using MS SQL Server Spatial database in PostgreSQL using foreign data wrapper?

texts: How to use MS SQL Server 2012 spatial database in postgresql using foreign data wrapper?


texts:
I generated this geojson data in my server:

{"type":"MultiLineString","coordinates":[[30.0212,40.123],[30.0229,40.1235],[30.0246,40.1239],[30.0262,40.1241],[30.027,40.1241],[30.0283,40.1241]]}

in ol3 I did:

// data above the client via ajax
$.each(data , function (k,v) {
     console.log(k,v);

      var g1=JSON.parse(v);

      var vectorFeature = new ol.Feature({
                geometry: g1 ,
      });

      featuresArr.push(vectorFeature);

})

var vectorSource = new ol.source.Vector();

vectorSource.addFeatures(featuresArr);

var yol2=new ol.layer.Vector({
            projection: 'EPSG:3857',
            source: [vectorSource]
        });

if( yol2 instanceof ol.layer.Vector) {
                  // alert("yol2");
                  olTrMap.addLayer( yol2);

               }

Then nothing is displayed.

Do you have any idea?


texts: Iterate through a list to use then use SearchCursor to populate new list

texts:
I have a featureclass of buildings.  Each one has a unique number ("BO_UNIQ").  I create a list of these unique numbers (1251 elements).  I then want to iterate though this list and get facility names from another table by building number. There are multiple facilities for each building number.  Eventually I will add the elements of the facility list to field(s) in the buildings featureclass.

As it is now the script will iterate one time with the da.SearchCursor

What am I missing?

# Import modules
import arcpy
from arcpy import env

# Set environment
env.workspace = r"PATH\EAD Facility Names\EAD_NAMES.gdb"
env.overwriteOutput = True

# set variables
buildList = []
fieldList = ["FACIL1", "FACIL2", "FACIL3", "FACIL4", "FACIL5", "FACIL6", "FACIL7", "FACIL8"]
eadBuild = "EAD_BUILDINGS"
eadTable = "tblEADFacilities"

# populate buildList with building numbers
with arcpy.da.SearchCursor(eadBuild, ["BO_UNIQ"]) as cursor:
    for row in cursor:
        buildList.append(row[0])
        print(row)

print(buildList)
print(len(buildList))

# Iterate through list of BO#s
for build in buildList:
    with arcpy.da.SearchCursor(eadTable, ["Name"], "BO =" + str(build)) as cursor:
        for row in cursor:
            facList = []
            facList.append(row)
            print(row)
    print(build)  # Will iterate and print all elements of the buildList and not run the SearchCursor branch.


texts: Does arcpy have a way to determine if a data frame overlaps a page extent?

texts:
I'm in the process of writing a script that will iterate through multiple MXDs and extract the data from layers that are turned on and only data that is within the data frame extent.

Some of the MXDs may have data frames that are not visible on the page. I don't want to extract the data from these data frames.  So, I want to write a function that will test whether a data frame is within the page extent.

My first thought was to access the dimensions of the page layout using the page size property of the Map Document. Then I would iterate through each data frame and test the X and Y position and compare them to the page dimensions.  In theory, this seemed like a good idea. However, I noticed that the positions change when the anchor point is changed.

Take the following example of an A4 page with four data frames.  Two on the page and two off the page:

Focusing on the data frame on the left, I can see that it's off the page.  When I look at the X position of the data frame, I can see that it is negative which seems to tell me that it is off the page.

Changing the anchor point, changes the X position:

and

I thought about adding the X position to the width.  In the first example, when the anchor is set to bottom left, the X position plus the width is still negative, which would tell me that the data frame is off the page.  However, if the anchor point is in the middle, or right, the sum of X position and width is positive, which would falsely indicate that the data frame is on the page.

This gets further complicated when the dataframe is partially overlapping the page. I also haven't gotten into the dataframe on the right or testing Y position, though the logic would be similar.

So, my question is, does there exist a method in arcpy to test if a data frame is on or overlaps a page extent?

I'm using ArcGIS Desktop 10.4.1 Basic.


texts: copy related table data into a single field in an attribute table

texts:
I have a featureclass of cities and a separate table of city administrators (multiple administrators for each city) and related contact info.

What I need to do is consolidate all city administrator info into the attribute table of the city featureclass. I'd like to add and populate fields for contact1, contactemail1, contact2, contactemail2, etc etc

A join, obviously, dumps all but the first record.

The relate shows me all the info i need, but I need to somehow copy all that info into a single row in the original featureclass, rather than clicking on each city to see the related table.

Possible?


texts: Losing raster values when exporting in QGIS

texts:
I have raster from Landsat 4 that I have created a symbology for. I want to save the new raster and to keep the NDVI values, but whenever I save the raster I lose the NDVI values.

On the left side you can see the values I get from the raster. I have saved band 4 which is supposed to be NDVI. But now it does not show the correct values anymore.


texts:
This page https://viewer.nationalmap.gov/services/ includes a link under "WCS Services" that lists a number of subdatasets.  This command
gdalnfo "https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer?request=GetCapabilities&service=WCS"

works fine.  I have been retrieving elevation data from the first subdataset:
SUBDATASET_1_NAME=WCS:https://elevation.nationalmap.gov:443/arcgis/services/3DEPElevation/ImageServer/WCSServer?version=2.0.1&coverage=DEP3Elevation
SUBDATASET_1_DESC=DEP3Elevation
and most of the time it works but I think I must be doing something wrong since sometimes there are large missing areas.
For instance this area (Boston) works just fine:
gdal_translate "WCS:https://elevation.nationalmap.gov:443/arcgis/services/3DEPElevation/ImageServer/WCSServer?version=2.0.1&coverage=DEP3Elevation" test.tif -projwin -71.07 42.36 -71.03 42.33 -projwin_srs EPSG:4326 -tr 100 100

but this command with only a different "projwin" (in Colorado) is about 75% missing:
gdal_translate "WCS:https://elevation.nationalmap.gov:443/arcgis/services/3DEPElevation/ImageServer/WCSServer?version=2.0.1&coverage=DEP3Elevation" test.tif -projwin -107.03 37.28 -107 37.25 -projwin_srs EPSG:4326 -tr 100 100

Can anyone help me understand what I'm doing wrong here and how to fix it?
I am using a very recent OSGeo4W featuring GDAL 3.0.4, released 2020/01/28.