texts: Using GDALwarp for reprojecting netCDF file?

texts:
I am trying to re-project a netCDF file using the gdalwarp command. However, the gdalwarp is giving me a black output.

Below is my code:

import gdal
import osr
import netCDF4
from osgeo.gdalconst import *
import os
import subprocess

gdal.AllRegister()
file_nc = 'input.nc'
ds = gdal.Open(file_nc, GA_ReadOnly)
Dataset = ds.GetSubDatasets()

substr = 'Rrs_655'
band = ''
names = [i[0] for i in Dataset]
for word in names[:]:
if word.endswith(substr):
band = word
b = 'gdal_translate -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 '+ band + ' output13.tif'
print(b)
subprocess.run(b, shell = True)
cmd = 'gdalwarp output13.tif outfile2.tif -t_srs "+proj=longlat +ellps=WGS84"'
subprocess.run(cmd, shell = True)

The gdalinfo of output13.tif band is:

Files: output13.tif
Size is 7821, 7931
Coordinate System is `'
GCP Projection =
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 7931.0)
Upper Right ( 7821.0, 0.0)
Lower Right ( 7821.0, 7931.0)
Center ( 3910.5, 3965.5)
Band 1 Block=7821x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Offset: 0.05, Scale:2e-06
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767

The gdalinfo if the file outfile2.tif is:

Files: outfile2.tif
Size is 12259, 12307
Coordinate System is:
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-20696.575935146218399,3354.169516545325678)
Pixel Size = (1.950162663734099,-1.950162663734099)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -20696.576, 3354.170) (Invalid angle,Invalid angle)
Lower Left ( -20696.576, -20646.482) (Invalid angle,Invalid angle)
Upper Right ( 3210.468, 3354.170) (Invalid angle,Invalid angle)
Lower Right ( 3210.468, -20646.482) (Invalid angle,Invalid angle)
Center ( -8743.054, -8646.156) (Invalid angle,Invalid angle)
Band 1 Block=12259x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767


texts: Joining Excel with ArcMap by ModelBuilder?

texts:
I am trying to use ModelBuilder in ArcMap 10.

I would like to join an Excel file (xls) with a shapefile.

The Excel file and the shapefile have one identical column with data but the Shp.file contains more data than the xls file.

Therefore I would like to make a join, based on the two identical columns, and get a shape file that contains the data that is identical?


texts:
I need to create a web map showing many raster layers. I am using mapserver and want to use Openlayers, however I can't find any good tutorial about it. I see a couple of old questions (1 and 2) saying that there were no good documentation. Have the things changed recently?

I would like to find a tutorial teaching from basics to rather complicated things with good explanations of the code and pictures/examples of the results. For now I managed to do only the simplest web map with my .map file, but I need to customize it (add legends, group layers, add more controls, embedding, etc).


texts: How to cast SDE table to IFeatureClassLoad

texts:
How do you cast an SDE table to IFeatureClassLoad using ArcObjects?  I'm getting an error with this attempt below which is not the correct method:

   Dim pSDEFeatWksp As IFeatureWorkspace = workspaceFactory.Open(propertySet, 0)
   Dim table2 As ITable = pSDEFeatWksp.OpenTable("database.LOADER.tablename")

   Dim fclass As IFeatureClass = table2 ' <<< ERROR HERE
   Dim FCLoad As IFeatureClassLoad = fclass
   FCLoad.LoadOnlyMode = False

ERROR

  System.InvalidCastException: Unable to cast COM object of type
  'System._ComObject to interface type
  'ESRI.ArcGIS.Geodatabase.IFeatureClass'...

I found this example which uses IFeatureClassLoad, however it is not dealing with a SDE.

This question is extending upon a previous question of mine.

Edit SDE table with ArcObjects - continued

Thank you!


texts:
Newbie in ArcObjects programming. Using C#. Developing a desktop application.

I need to join 2 tables:

POI layer which is a Point layer and
a parcel layer which is a Polygon

Where

POI.HouseNo = Parcel.UNIT_ID

and then create a new Shapefile which has only the 4 fields:

POIName,
POIType,
HouseNo (all from POI layer attributes),
ParcelName (all from Parcel Layer attributes).

The user selects the POI Layer and the Parcel layer from combo boxes.

Here is the code I have written so far:

IFeatureWorkspace featureWorkspace;
IQueryDef queryDef = featureWorkspace.CreateQueryDef();
queryDef.Tables = cboPOI.Text + "," + cboParcel.Text;
queryDef.SubFields = cboPOI.Text + ".Name," + cboPOI.Text + ".Type," + cboPOI.Text +".HouseNo," + cboParcel.Text + ".ParcelName";
queryDef.WhereClause = cboPOI.Text + ".HouseNo =" + cboParcel.Text+ ".UNIT_ID";

Hope I'm on the right track.
Now I need to run this query and get the output to a new Shapefile. How can I do this?

please help.


texts: Calculating Azimuth between two points in PostGIS?

texts:
I have a GIS API based on OpenLayers. I have been trying to implement Azimuth calculation in JavaScript and I needed some way to calculate the Azimuth in order to perform tests.

I started using PostGIS, but there seem to be many ways to calculate the Azimuth between two points. I show you three of them and, some return different results.

-- Example 1 - Result 90

select st_azimuth(
st_transform(st_geomfromtext('POINT(-81328.998084106 7474929.8690234)', 900913), 4326),
st_transform(st_geomfromtext('POINT(4125765.0381464 7474929.8690234)', 900913), 4326)
)/pi()*180

-- Example 2 - Result 155.692090425822
SELECT degrees(
ST_Azimuth(
ST_MakePoint(-81328.998084106, 7474929.8690234)::geography,
ST_MakePoint(4125765.0381464, 7474929.8690234)::geography)
)

-- Example 3 - Result 90

SELECT degrees(
ST_Azimuth(
ST_MakePoint(-81328.998084106, 7474929.8690234),
ST_MakePoint(4125765.0381464, 7474929.8690234))
)

What is the correct way to calculate the Azimuth in PostGIS?

Of course, I want to consider geodesic coordinates.

Is there a way to know how PostGIS realizes the calculations?

I mean, is it possible to see the way the "ST_Azimuth" function is implemented?


texts:
What is the easiest way to deliver a map completely offline? E.g. like if I wanted to package it up and put it on a mobile device that is not in range of wifi.

How would one serve up map tiles w/o a tiling server?

Specifically looking for something javascript based.


texts: Calculating area of raster pixel when pixels are different sizes at different latitudes using ArcGIS Desktop?

texts:
I have a large raster file of total mangrove cover for the whole globe (obtained from NASA: http://sedac.ciesin.columbia.edu/data/set/lulc-global-mangrove-forests-distribution-2000/metadata). I'm trying to calculate the area of mangrove forest that fall within protected areas (that are defined by polygon layer).

I can calculate the number of pixels that fall within the polygon, however, as the area of the pixels changes at differing latitudes, how can I reliably convert these pixel counts into area?

I'm aware that flattening a sphere causes some sort of distortion, so I tried Define Projection tool to change to World Cylindrical Equal Area (from WGS84), but that just shrinks the whole raster down to an unusable scale.

If I set the coordinate system of the layer to World Cylindrical Equal Area before importing the raster, it preserves the scale, but after measuring the pixels with a ruler, I can see they change at different latitudes.

I've included some screenshots for clarification.

The first image is the global mangrove forest cover (dark green), the second and third images are the >Layer Properties >Source window of the mangrove.tiff raster file (unaltered). And the last image is a measurement of 1 pixel close to the equator.


texts: Problems setting up eclipse (with PyDev) for ArcGIS 10

texts:
I'm trying to use eclipse with PyDev to program in python. In Windows 7 x64, I've added a PATH environmental variable (PATH C:\Python26\ArcGIS10.0\python.exe;) and what I thought were the relevant folders to PYTHONPATH (C:\Dropbox\CODE\Python\; C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy)

However, when I type the command:

import arcpy

in eclipse I get an error message. However, when I try importing ArcPy via IDLE I have no problems. Is there some way that I should be including arcpy in the build path that I am missing?


texts: Union Tool returns an not-before-seen error in QGIS?

texts:
'NoneType' object has no attribute 'isGeosEmpty'

Above the error I get when attempting to combing two polygons of the same spatial extent, but of course with two different sets of information. Both have same projection defined. not sure why this isn't working. Do the attribute tables have to have at least one column of similar information?


texts: Adding temporary graphic or vector layer to QGIS to replace the slow select feature method

texts:
I have a custom QGIS 1.9 dev (so 2.0 I guess) plugin which collects and displays a list of polygons in a QTreeWidget. I would like to select a list item and it immediately display a highlighted polygon to which it corresponds. This is possible by using the QgsVectorLayer.select() method, but the method is very slow and I would like it to display almost instantly. I can simulate the highlighting of a polygon by creating a new in-memory vector layer and adding a feature using the geometry of the selected list item, but I don't want to display that new layer in the Layers window.

My code (not 100% usable here) shows the basic implementation so far. My preferred outcome would be a fast 'select' method, or else a way to display a new temporary vector layer without displaying it in the Layer window.

class FindLandParcelDialog(QtGui.QDialog):
    def __init__(self, iface):
        QtGui.QDialog.__init__(self)
        self.ui = Ui_FindLandParcel()
        self.ui.setupUi(self)
        self.parcel_tw = self.ui.treeWidget

        self.parcel_dict = {CP100: (143, 'shape_wkt'), CP101: (144, 'shape_wkt')}

        self.parcel_tw.itemSelectionChanged.connect(self.selectParcel)  # Very slow select
        # self.parcel_tw.itemSelectionChanged.connect(self.highlightParcel)  # displays new layer in Layer window - not wanted

    def selectParcel(self):
        '''Selects polygons with the ID stored in the dictionary linked to the
        QTreeWidget item'''
        self.cLayer = self.canvas.currentLayer()
        idx = self.getTreeCurrentItem(self.parcel_tw)
        parcel_id = self.parcel_dict.items()[idx][1][0]

        self.cLayer.removeSelection()
        for feat in self.cLayer.getFeatures():
            if feat.attributes()[0] == parcel_id:
                self.cLayer.select(feat.id())

    def getTreeCurrentItem(self, QTreeWidget):
        '''Returns index of the selected top level QTreeWidgetItem'''
        currentitem = self.parcel_tw.currentItem()
        idx = QTreeWidget.indexOfTopLevelItem(currentitem)
        return idx

    def highlightParcel(self):
        self.cLayer = self.canvas.currentLayer()
        self.zoomtolayer_btn.setEnabled(True)
        idx = self.getTreeCurrentItem(self.parcel_tw)
        parcel_wkt = self.parcel_dict.items()[idx][1][-1]
        self.createInMemoryFeature(parcel_wkt)

    def createInMemoryFeature(self, wkt):
        layer = QgsVectorLayer('Polygon?', 'temporary_polygons', 'memory')
        pr = layer.dataProvider()
        layer.startEditing()
        feat = QgsFeature()
        feat.setGeometry(QgsGeometry.fromWkt(wkt))
        pr.addFeatures([feat])
        layer.commitChanges()

SOLUTION The working code to display a temporary feature (from Nathan W's answer):

class FindLandParcelDialog(QtGui.QDialog):
    def __init__(self, iface):
        QtGui.QDialog.__init__(self)
        self.ui = Ui_FindLandParcel()
        self.ui.setupUi(self)
        self.parcel_tw = self.ui.treeWidget

        self.parcel_dict = {CP100: (143, 'shape_wkt'), CP101: (144, 'shape_wkt')}

        self.parcel_tw.itemSelectionChanged.connect(self.highlightParcel)

    def highlightParcel(self):
        self.cLayer = self.canvas.currentLayer()
        idx = self.getTreeCurrentItem(self.parcel_tw)
        parcel_wkt = self.parcel_dict.items()[idx][1][-1]
        self.createRubberBandFeature(parcel_wkt)

    def createRubberBandFeature(self, wkt):
        band = QgsRubberBand(self.canvas)
        band.setToGeometry(QgsGeometry.fromWkt(wkt), self.cLayer)
        band.setColor(QtGui.QColor(255,255,0,255))
        band.show()

    def getTreeCurrentItem(self, QTreeWidget):
        '''Returns index of the selected top level QTreeWidgetItem'''
        currentitem = self.parcel_tw.currentItem()
        idx = QTreeWidget.indexOfTopLevelItem(currentitem)
        return idx


texts: Google Earth Engine: change multiple labels in legend for time series by region graph

texts:
I am producing a time series graph comparing the mean NDVI in two different types of regions. Both types exist out of multiple polygons (MultiPolygons). However, I keep failing to name these two regions in my legend of my time series graph. Adding it as a 7th argument in ui.Chart.image.seriesByRegion only works if I have one region since it only accept one string and not a list of strings. How can I change labels in the legend of a graph of series (as indicated .setSeriesNames did not work for me either).

// Create a time series chart.
var regions = ee.FeatureCollection([Small_Avocado_Farms,Large_Avocado_Companies]) // featurecollection from two multipolygons
var tempTimeSeries = ui.Chart.image.seriesByRegion(
    all_images_ndvi, regions, ee.Reducer.mean(),'NDVI',30)
        .setChartType('ScatterChart')
        .setSeriesNames(['Small','Large'])
        .setOptions({
          title: 'NDVI of small- vs. large-scale avocado farms ',
          vAxis: {title: 'NDVI'},
          pointSize: 4,
          legend: {position: 'top'}
          })

print(tempTimeSeries)

The following answers did not work for me:
Change the label of the ui.Chart.image.seriesByRegion in Google Earth Engine


texts:
I would like to change the splash screen of QGIS 3.4 LTR. I read the instructions from Customise the splash image of QGIS but I didn't manage to change it.

Does the overwrite of the splash screen also work with QGIS 3.4? I tried to add:

[Customization] splashpath=C:\Temp\QGIS-3.4.8-1-64bit\splash.png

to QGISCUSTOMIZATION3.ini, but the splash screen is replaced by with square.

I guess it must be an issue related to the path, must it be relative? Or is the size of the picture important?


texts: Using GDAL in Python to Convert NetCDF to ASCII and flip latitudes

texts:
I have a NetCDF file containing bathymetry. I am using the GDAL and NETCDF4 libraries to export the file to ESRI ASCII using the 'AAIGrid' format. However, a common problem is the latitudes are in the reverse order (see here for patch made to gdal and here for Matlab and here for R solutions).

I think the solution can be achieved using Numpy and flipud function. I am not sure whether to flip the bathymetry or the latitude. When plotting using matplotlib, both seem to work to fix the problem.

I seem to be encountering two problems. First, when using gdal.Translate, the output Dataset cannot be passed as the NetCDF4 'Dataset' object - it must be a file. I would like to pass it a Dataset object which I have modified by flipping the latitudes. When I pass in an object (even an unmodified Dataset object which just references the original .nc file), I get the following error:

#gdal.Translate("out.txt", "test4.nc", format="AAIGrid", noData="-10")
gdal.Translate("out.txt", new_netcdf, format="AAIGrid", noData="-10")

    TypeError: in method 'TranslateInternal', argument 2 of type 'GDALDatasetShadow *'

Second, I cannot seem to flip the latitudes by modifying the original file.

from netCDF4 import Dataset
from osgeo import gdal
import numpy as np

file = "ORIGINAL\Bathymetry_1km.nc"
netcdf_f = Dataset(file, "r", format="NETCDF4")

latitudes = netcdf_f.variables['lat']
lat_flipped = np.flipud(latitudes)
netcdf_f.variables['lat'][:] = lat_flipped

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-126-457cda137417> in <module>
----> 1 netcdf_f.variables['lat'][:] = lat_flipped
netCDF4\_netCDF4.pyx in netCDF4._netCDF4.Variable.__setitem__()
netCDF4\_netCDF4.pyx in netCDF4._netCDF4.Variable._put()
netCDF4\_netCDF4.pyx in netCDF4._netCDF4._ensure_nc_success()
RuntimeError: NetCDF: HDF error

So, next I create a new NetCDF 'Dataset' object (actually as a file on the drive because I can't figure out an alternative) and populate it with the lats, lons, and bathy dimensions / variables from the original dataset. That seems to work except when I go to translate it to ASCII, I end up with a grid full of 9.9692099683868690468e+36 and nothing else.


texts: Geoserver KML Rendering in Google Earth - SLD Errors

texts:
I have a layergroup with custom SLD's in Geoserver:
Layer preview

All the symbology renders fine. If I try and add the KML as a network link into google earth using this URL I lose most of my SVG and Polygon symbology.

The svg files are in the styles folder in geoserver and I've also tried copying one of them to an accessible URL. The SLD for that layer now contains:

<PointSymbolizer>
        <Graphic>
          <ExternalGraphic>
            <OnlineResource xlink:type="simple" xlink:href="http://www.dassh.ac.uk/images/carpark.svg" />
            <Format>image/svg+xml</Format>
          </ExternalGraphic>
          <Size>
            <!--<ogc:Literal>12</ogc:Literal>-->
            12
          </Size>
        </Graphic>
      </PointSymbolizer>

And as before, displays fine in OpenLayers (within Geoserver) but not Google Earth.

Any help gratefully received.

Thanks

Dan


texts: Finding all segments that makeup a road in OSM

texts:
I've spent 2 days trying to figure this out and I'm afraid that it's something that's not possible.  If that's the case then I'm thinking the wrong with with respect to OSM data.

Let's say I have a street, called Main Street.  On the displayed map on openstreetmaps.org it looks like one road.  To the underlying data it is made up of multiple ways (and is part of a relation of a larger route.)

Let's say I am only interested in the segment that passes through my town which is still made up of multiple ways.

Is there a (easy) way to find out what the ways are which make up this segment?

Is there a way to figure out ways that touch (either the beginning or end) of a way?

I can think of an expensive way (in terms of resources) doing this but I'm afraid that I'm looking at this the wrong way.

I was initially using a XML file that I downloaded of my town but I've also loaded the data into a postgre database (with postgis).

I've read a bunch of questions on stackexchange and the osm website as well.
I've looked at this as well: http://wiki.openstreetmap.org/wiki/Elements


texts: Need help regarding Hot-spot Analysis through ArcGIS JSAPI 3.6/3.7

texts:
I am going through Hotspot Analysis through java script sample and I want to perform similar thing at my end.

While accessing this sample it taking too much time to display the result (for sign in I created 30 day trial account so is that matter to access the sample?)

-According to this sample ESRI mention that if you want to use this type of service then you need to subscribe the ArcGIS Online analysis service (Why ? we have ArcGIS  Server which itself costly product then why should we again purchase the ArcGIS Online?)

-Can I create similar service at my end using geo-processing and use in the code?

-In future if I want to use spatial analysis services then ArcGIS Online subscription will be necessary ? (If so then whats the importance of having ArcGIS Enterprise Advance license?)

Any help regarding this will be great !

Software suit: ESRI 10.1 and ArcGIS JSPAI version 3.6/3.7


texts:
In ModelBuilder I am trying to perform the following query, using the select by attribute tool. I am having great difficulty with the sql syntax.

What I am trying to do is

Select Top 3 (Avg) Price, GasID Group by GASID order by GASID

The help indicates that aggregrate functions need to be performed in subqueries only?
Is this true.

I don't want a subquery.
Only this works but returns no rows and it is not what I am looking for:

"RegGasPrice" = (SELECT MAX( "RegGasPrice") FROM gasprices).  There is no "in" so I must use equal.


texts: Using Mapproxy parameters as input to Mapnik data source

texts:
I'm using a structure of Mapproxy and Mapnik to render raster data through a WMS.
I was able to connect the Mapnik source on Mapproxy and render a unique file from my XML Mapnik data source configuration.
I would need to specify a URL parameter to specify the data the user requested. So, I would like to create a dynamic script that can receive parameters from the URL and process the raster output according to those parameters.
Is there a way to set up the Mapproxy to invoke a Mapnik Python script or an API, instead of invoking the Mapnik XML?


texts: Filling hole in DEM using ArcGIS Desktop?

texts:
How would I fix my DEM because it has holes?

I need to maybe extrapolate from surrounding terrain.

Here is the picture


texts:
I added a new column in my PostgreSQL/PostGIS database for getting surface area. I can open as a layer in QGIS, but when I draw a new polygon the area column is not updated, is there a way to have it updated as soon as I finish drawing the polygon?

This is the SQL I used:

ALTER TABLE public."EDIFICIO_2" ADD COLUMN area double precision;
UPDATE "EDIFICIO_2" SET area = round((ST_AREA(geom)) :: numeric,2);


texts: represent multiple geometry as single feature

texts:
There is a school. it is on two disjoint areas, one is primary wing and other is secondary. I need to make single feature from the two areas such that, there is only one column in gis table, i.e. only one record.

Also is it possible for them to have separate values but still show up as one record.


texts: Is there a way to set the Geocortex Essentials in full-screen mode by default?

texts: I was wondering if there is a setting in Essentials that would allow me to publish my webgis portal in full-screen mode by default?


texts: Unable to render layers in MapCanvas using PyQGIS/PyQt4

texts:
I am trying to get some layers into a MapCanvas and style them.  I am following the PyQgis Cookbook info here but my code is just not working.  I've been spinning my wheels on this for over a week now and there's no one in my office familiar enough with python/qgis.  What am I missing???

QgsMapLayerRegistry.instance().addMapLayer(layer1)
canvas1.setExtent(layer1.extent())
canvas1.setLayerSet( [QgsMapCanvasLayer(layer1)] )
Traceback (most recent call last):
File "<input>", line 1, in <module>
NameError: name 'QgsMapCanvasLayer' is not defined

In R you just do this:

library(rgdal)
nta <- readOGR("my/path", "mylayer")
spplot(nta)
library(sp)
spplot(nta, zcol="symbolAttributes")


texts:
I have LAS files with code 0 - never classified. I have no other classes other than this.

How do I assign class codes for my LAS file?


texts: proj4: how to convert UTM to local coordinate system

texts:
I'm trying to convert a set of UTM32N coordinates to a local coordinate system with the same datum.

I have the four Helmert parameters (x,y,scale,theta) that defines the coordinate system.

How can I apply this transformation using proj4 strings?

I've tried the following one but it differs from the expected result:

proj +proj=helmert +ellps=WGS84 +x=-711066.2010726777 +y=-4532044.137231419 +s=0.999785911377552 +theta=-2.089700973339795


texts: Unable to use GeoJSON using leaflet package in R

texts:
I am trying to plot the GeoJSON in R using Leaflet package. Below is the code and the error.

library(geojsonio)
library(leaflet)
library(data.table)
library(plyr)
library(rgdal)
library(sp)
library(RColorBrewer)
library(rgeos) #for simplification
library(leafletR)
library(sf)
library(jsonlite)
library(RJSONIO)

mydata <- fromJSON("https://gist.githubusercontent.com/senthilthyagarajan/eb7a2771eab4639e94d5f9eaad28cb33/raw/1cfe355a56d2c1856a70a5389a4eadf06d782748/data.geojson",flatten=TRUE)

leaflet(mydata) %>%
  addPolygons(color = "#444444", weight = 1, smoothFactor = 0.5,
              opacity = 1.0, fillOpacity = 0.5,
              fillColor = ~colorQuantile("YlOrRd", nghbrhd)(nghbrhd),
              highlightOptions = highlightOptions(color = "white", weight = 2,
                                                  bringToFront = TRUE))

Error: lexical error: invalid char in json text.
                                       FeatureCollection
                     (right here) ------^

Please ignore the long list of packages mentioned above.


texts: How do I get revised data from the U.S. Census?

texts:
I need to calculate population density at the level of U.S. Census tracts for across a number of states.

I used the instructions here Where to get 2010 Census Block data? to download DP-1 for the states CT, ME, MA, NH, NY, RI, and VT, but some of the values in the column for Total Population show up like this:

2114(r39570)

This evidently means that the value has been revised.  I went back to American FactFinder to view the data, and by clicking on the value "(r39570)", I'm able to confirm that these indicate that the value has been revised:

As far as I can tell, the correct value, which is evidently 1120 in this example, does not appear anywhere in the .csv that I imported.

How can I get the correct values onto my computer and appropriately imported into QGIS so that I can perform the correct spatial calculations?


texts: How do I activate "Open GRASS tools" button?

texts: I have opened two GRASS layers and need to use the "Open GRASS tools" button. This button remains unactivated somehow, while the "Edit GRASS vector layer" right next to it works perfectly fine. What can I do?


texts: unable to find "attribute" field in "ClassBreaksRenderer " class

texts:
I am trying to run Map ServerLayer sample given on arcgis resource center.When I pasted the code in my FlashBuilder 4.6 then it is showing error on this line <esri:ClassBreaksRenderer attribute="POP2000"> . Error is

Multiple markers at this line:
-Cannot resolve attribute 'attribute' for component type com.esri.ags.renderers.ClassBreaksRenderer.
-ClassBreaksRenderer
-1 changed line

It means it is not recognizing the attribute field of esri:ClassBreaksRenderer. I am unable to understand the reason?Any suggestions?

Note: I am using flex4.6.0 and latest version of arcgis API which is agslib-3.3-2013-04-24.swc


texts: Reconcile/POST version DB edits- only by selected USERS(last edited)

texts:
I have been reading a bunch of helpdocs/whitepapers from ESRI and could not find a solution.
My Scenario:
I have a multi-user DB and want child 1 edits concile/POST to parent 1 which then parent 1 All this works fine however I foresee a potential issue.

Child 1 version is for field crews making edits/creating data in the field.
parent1 is for our team manager to Q/A these edits which they will then reconcile/post to their parent 1 version. From there I have a script that runs nightly to reconcile/post parent1 to default.

I have created a simple script tool for the manager to execute once they verify all the edits which will only post child1 edits to parent1.

However what if there were 50 edits made and the manager only wants to push 20 of those edits to parent1 (either by username or by individual OID(row) by table). Is there a way to select only certain rows to reconcile/post perhaps by last user to edit (editor tracking is enabled) which I can implement into my script as a GetParameterasText().

I basically need to modify the reconcile manager arcpy command by only selected users or by OID by table/layer.  From what I have read it seems to be a all or nothing for reconcile/POST
Hope this makes since if not I will try to re-write question.


texts: MapInfo 11.5 - Legend Designer - problem with symbols

texts:
First of all I must say I love new legend designer. It is great to have finally better tool for positioning and even some changing the size of lines and region symbols, hooray. Big step for MapInfo to be one day user friendly professional tool for map layouts.

But... when I use Legend Designer all content of legend window is automatically converted to bitmap for layout print/export to PDF. Is there some option to turn this off and keep text/vectors? I found out this way I can't use it at all with symbols.


texts:
I have a DEM of Saint Lucia, obtained from the CHARIM database. It is a large geotiff (150MB), and I would like to crop it as I only need a small part of it. I have been trying to use Raster -> Extraction -> Clipper in QGIS, but it fails because it says I need thousands of TB of disk space to complete the operation. From reading around, I presume that this is because the file is already compressed, however I don't really know how to go about solving it.

This is the layer metadata:

GDAL provider
GTiff
GeoTIFF
Dataset Description
/Users/cameron/Desktop/Design IV/St Lucia GIS/topography.tiff
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=1 (unitless)
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Band 1
STATISTICS_MAXIMUM=779.13000488281
STATISTICS_MEAN=nan
STATISTICS_MINIMUM=-inf
STATISTICS_STDDEV=nan
Dimensions
X: 4536 Y: 8763 Bands: 1
Origin
-61.0813,14.1115
Pixel Size
4.61889e-05,-4.61889e-05
No Data Value
-1e+38
Data Type
Float32 - Thirty two bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
-61.0812601912316993,13.7067470580471564 : -60.8717473436720837,14.1115003832598784
Band
Band 1
Band No
1
No Stats
No stats collected yet

Can anyone tell me what I should be doing instead? I am relatively new to GIS and this is the first thing I have come across that I can't immediately find my way around.

EDIT: The Clipper settings are as follows:

gdal_translate -projwin 507118.619912 1549423.97735 511757.851837 1546542.77015 -of GTiff

And the error log:

Warning 1: Computed -srcwin 1.09806e+10 -3.35451e+10 1.0044e+08 6.23788e+07 falls completely outside raster extent. Going on however.
ERROR 3: Free disk space available is 8291942400 bytes, whereas 25061401106245716 are at least necessary.`


texts: Loop to Apply Multiple Symbology to Multiple Layers

texts:
I am trying to apply multiple symbology to multiply layers that I have created from multi-tasking custom tool. The tool applies buffer, clip tool to multiple feature classes and then all the feature classes are saved as feature layer.

I have also pre-defined .lyr files. This involves using  arcpy.ApplySymbologyFromLayer_management(in_layer,in_symbology_layer)

I would like to apply these pre-defined .lyr files to the feature layers that have been created. Pre-defined and created feature layers have the same name.

Example:

This code I have is as follows:

import os
import arcpy

filelist = os.listdir (r"C:\gdb\Created_Layers")

layerlist = os.listdir (r"C:\Tool\PreDefined_Layers")

for f in filelist:
    for s in layerlist:
       if f == s:
        arcpy.ApplySymbologyFromLayer_management(f,s)

But it doesn't work. I still see the random defined symbology for every feature layers that have been created.


texts: How to label a dynamic KML file in OpenLayers?

texts:
How can I display labels for features in a KML file on an OpenLayers map. I was thinking of using popup messages, but labels are more attractive and eye-catching.

I want to replicate the labeling method of Google Earth:

Any hints? Thanks in advance.


texts:
I have 11 Satellite images in .jpg format showing the same area but in 11 different years. They align perfectly and have the exact same extent. They are not georeferenced and have 3 bands each.
Is there a way to georeference them all at the same time?
My idea was to 'stack' them first, georeference them and 'unstack' them again.
That is how I want to avoid georeferencing the images one by one.
I made a VRT with QGIS but I don't know how to undo that VRT later so that I have the single rasters again. With ArcGIS I tried to mosaic them but it seems that it is not the right way.
Any idea?


texts: Moving nodes of different multipolygon on different layers

texts:
I would like to import into Openstreetmap the land-use data. In order to do a good job, I would like to correct some imperfections before the import. Take a look at this screenshot (which is a small area of the current project):

The two green regions represent two different adjoining towns. Each multipolygon indicates a land use. The red lines are main roads, tracks, paths, water streams. It is simple to recognize main roads, since these separates the land-use with a white space.

Unfortunately, in the picture you can see the main road is going to cross the separation regions on multiple occasions. The main roads are accurately mapped with GPS devices, while the land-use data has been built by local administrations knowledge and aerial imaginery: this explains the crossing.

I would really like to solve all these crossing before the import. Given my limited experience with QGIS, the only solution I can see is to manually move all the different nodes of the different multipolygon and the different layers (with the help of the snap options): this is simply unacceptable given the tremendous amount of time it would take.

Do you know any other possible approach? Is there any plugin that could help me achive my goal?


texts: Verifying that attribute indexes are being used in query from ArcObjects?

texts:
We have some software as part of a larger project written using ArcObjects via .NET. This software migrates and merges data from multiple file geodatabases in one schema to a single file geodatabase in a similar but different schema. Some steps in the migration involve looking up an existing row and updating it with data from a different table. This can get very slow when the destination geodatabase is large, so I have added some multi-field attribute indexes to improve the query/update performance.

Is there a way that I can confirm that the queries are using the new attribute indexes, and perhaps how the query is being executed?

I guess I'm really looking for the equivalent of a query plan - something similar to EXPLAIN PLAN in Oracle.

The queries are on a file geodatabase using ArcObjects directly (no RDBMS backend or SQL).


texts: Tutorial on creating geospatial visualizations with Processing(.org)?

texts:
Recently, I've seen examples, (lots of examples ) which use Processing(.org) for creating network visualizations. So far though, I have not been able to find any good tutorials on it.

I was wondering if anyone knows a good resource for learning this specific type of visualization?


texts: Shapefile to Rectangle Overlay/Intersect Test?

texts:
I have a shapefile and a repository of images that have corner Lat/Long points. I want to search that images who overlaps/intersects to shapefile.

Here is the scenario that I am doing, creating a polygon from each image corner lat/long points and placing it in temporary memory.

Now how can I check whether this created polygon intersects/overlaps the shapefile?

Any solution recommended for this issue?


texts: arcpy perform function on features in feature class

texts:
My data consists of a directory of shp files.
I'm trying to accomplish the following:
   Go through each shp file and perform a function on each feature.
The way my script is working now, it is preforming the function on every feature in the feature class at once, so the fields I'm writing to has the field values of every feature in the feature class in it rather than the values of each feature seperately. Here is my script as of now:

import arcpy, time, datetime, os
from datetime import datetime
from arcpy import env
arcpy.env.workspace = r"C:\Data_State\MT\temp\gpx"
def PopulateGPX():
    GPXFile = r"C:\Data_State\MT\temp\gpx\GPXtest.gpx"
    writeGPX = open(GPXFile, 'a+')
    writeGPX.writelines('''
    <wpt lat="''' + str(Latitude) + '''" lon="''' + str(Longitude) + '''">
        <name>''' + str(Name) + '''</name>
        <cmt></cmt>
        <extensions>
            <gpxx:WaypointExtension>
                <gpxx:Address>
                    <gpxx:StreetAddress>''' + str(Address) + '''</gpxx:StreetAddress>
                    <gpxx:City>''' + str(City) + '''</gpxx:City>
                    <gpxx:State>''' + str(State) + '''</gpxx:State>
                    <gpxx:PostalCode></gpxx:PostalCode>
                </gpxx:Address>
            </gpxx:WaypointExtension>
        </gpxx:extensions>
    </wpt>
    '''+ '\n')
    writeGPX.flush()
for fc in arcpy.ListFeatureClasses():
    for feature in fc:
        Latitude = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "Latitude")]))
        Longitude = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "Longitude")]))
        Name = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "Name")]))

        Address = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "ADDRESS_TA")]))
        City = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "CITY_TAX")]))
        State = sorted(set([r[0] for r in arcpy.da.SearchCursor(fc, "STATE_TAX")]))

        PopulateGPX()

I hope I explained what I'm trying to accomplish well enough. Any ideas of where I'm going wrong? Using 10.3


texts: Different angle of the same Data set in different map layouts

texts:
I am trying to create a new layout - to get the bigger area of the map.
I used the same data: raster, vector.
Even if north arrow showing on the north in both cases:
there is a difference:

The original:

new one:

What can be a problem?


texts: Create a loop to run and store the r.stats output of several viewsheds?

texts:
I have 15 viewsheds (view 0-14) from which I would like to export the cell count. I can do this for one viewshed through:

count=[]
i=0

viewshedSize=grass.read_command("r.stats", overwrite=True,flags="c", input="view0")
countsperline= viewshedSize.split("\n")
count.append(type(float(countsperline[1].split(" ")[1])))

How can I run a loop to do this for each viewshed?


texts: ArcGIS layer to KML file contains  how can I get KML file without it?

texts:
So when I create line dot ect. after converting layer to KML file contains  "MultiGeometry" tag. How could I get KML file without it?
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
 xsi:schemaLocation="http://www.opengis.net/kml/2.2 http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd http://www.google.com/kml/ext/2.2 http://code.google.com/apis/kml/schema/kml22gx.xsd">
<Document id="Untitled_map">
  <name>Untitled map</name>
  <Snippet></Snippet>
  <Style id="PolyStyle00">
    <LabelStyle>
      <color>00000000</color>
      <scale>0.000000</scale>
    </LabelStyle>
    <LineStyle>
      <color>ffa5e1c5</color>
      <width>5.925750</width>
    </LineStyle>
    <PolyStyle>
      <color>ffa5e1c5</color>
      <outline>1</outline>
    </PolyStyle>
  </Style>
  <Folder id="FeatureLayer0">
    <name>Untitled layer</name>
    <Snippet></Snippet>
    <Placemark id="ID_00000">
      <name>palogonas</name>
      <description><![CDATA[]]></description>
      <styleUrl>#PolyStyle00</styleUrl>
      **<MultiGeometry>**
        <Polygon>
          <extrude>0</extrude>          <altitudeMode>clampToGround</altitudeMode>
          <outerBoundaryIs><LinearRing><coordinates> SOME COORDINATES </coordinates></LinearRing></outerBoundaryIs>
        </Polygon>
      **</MultiGeometry>**
    </Placemark>
  </Folder>
</Document>
</kml>


texts: Getting tiles from own Tilestache server to Openlayers application

texts:
The title sums up my current problem  -

I have installed and set Tilestache in our webserver and the tiles are accessible in browser
My developer is working on an application to overlay some custom layers in  kml format , over some basemap. We were using a basemap from Openlayers example, but now want to replace the url with tiles from the Tilestache.

Which method we have to get the tiles ? XYZ ? UTF grids ? httprequest ? Experts please guide..


texts: Formatting scale bar text with thousands separator in QGIS

texts: I wonder if this is possible. I want to show the text (numbers) of the scale bar in Print Composer with a thousand separator in order to make the labels more readable. All I can do for now is change the font and color. I can not find an option similar to the grid labels, where you can define a custom setting for the label and use format_number(@grid_number,0 for instance. I am using QGIS 2.18 in Win 10.


texts:
I am getting a pretty interesting error in my attempt to import arcpy in my python command line. The end goal of this effort is to have a working Spyder IDE that can run arcpy. I've been able to get this to work with the 10.5 version of ArcGIS by setting the correct paths but am unable to replicate it. Here's a rough outline of the steps I've followed:

32 bit Python2.7 downloaded
Anaconda Navigator downloaded and created an environment for Arc
Set Desktop10.6 path from 'C:\Python27\ArcGIS10.5\Lib\site-packages' to 'C:\Anaconda2\envs\arcgis10-5-1\Lib\site-packages.'
Created path 'C:\Python27\ArcGIS10.5\Lib\site-packages' in text editor and saved as 'C:\Anaconda2\envs\arcgis10-5-1\Lib\site-packages\
arcgis-python-site-packages.pth.'

Then I simply went to the Python command line to 'import arcpy' and here is my error:


texts: How to Convert from Point to Continuous Raster in ArcGIS 10.1?

texts:
I have a vector point dataset that I would like to convert to raster. The point-to-raster conversion tool creates a raster grid with spaces between each grid cell. I want to create a 'continuous' raster whereby each grid cell expands outwards until its' borders meet with the adjacent cell's borders.

I find the interpolation tools under spatial analysis to be unsatisfactory (spline is the closest I come to satisfactory). I do not want to smooth my data when I put it in raster form (ie I want to keep peaks and valleys in data).

Is there a geoprocessing tool that gets at what I am attempting to do?


texts:
I'm trying to convert this TIF to EPSG:3857, but when I run:

gdalwarp -of GTiff -s_srs epsg:4326 -t_srs epsg:3857 -co "tfw=yes" original.tif converted.tif

I just get a black TIF do you know what I am doing wrong?

Command output

> gdalwarp -of GTiff -s_srs epsg:4326 -t_srs epsg:3857 -co "tfw=yes" original.tif converted.tif
Creating output file that is 4088P x 4105L.
Processing input file original.tif.
Using band 2 of source image as alpha.
0...10...20...30...40...50...60...70...80...90...100 - done.


texts: Using MMQGIS to create accurate buffers?

texts:
I have read in several places that the best way to create a buffer around a vector point is using the MMQGIS plugin.

I want to create buffers in London locations within a few metres of each other which meant that Vector > Geoprocessing tools > Buffer is not accurate enough using decimal degrees.

I've downloaded MMQGIS and I understand how it works in principle, however when I create my buffer, even of only 1km it ends up sitting as one big lump west of South Africa and not in London.

I have checked the CRS for each of my layers which is WGS84 (EPSG:3857). I need to use this as I'm working with a unique historic map which has been formatted to this CRS.

When I checked my project CRS (Project > Properties > CRS) the project is set to Generated CRS (User:100000). I changed the CRS to WGS 84, but this just warped the map and the buffer was still west of South Africa. I have therefore changed the project CRS back to Generated CRS.

Could anyone explain to me the MMGIS fields and where I may be going wrong?

Or an alternative programme I could use to get my buffers down to metre accuracy?

I'm quite IT literate, but coding etc is a little beyond me.

I am working with QGIS3.2.3-Bonn on Mac.


texts: Representing degrees between -1 and 0 in Degrees:Minutes:Seconds

texts:
When representing latitude and longitude in degrees, minutes and seconds, with

-180 <= D <= 180 (or -90 and 90)
0 <= M < 60 (0 if D is at an extreme)
0 <= S < 60 (0 if D is at an extreme)

How can values between -1 degree and 0 degrees be represented? Since there is no distinction between 0 and -0, there is no way to distinguish between 0:12:34 and -0:12:34.

Is this just a limitation of this kind of representation, or is there some convention that allows it to cover this gap?


texts:
This is a WMS layer i'm working with:

It shows a list of various line data, 402-914.  As you can see to the right, this corresponds with various info on a topographic map.  My problem is that all I want are the contour lines, not everything else.  Unfortunately, I can't seem to find a way to separate all the different schema.  This is what the WMS page looks like:

They break it out by default into 4 different packages, each with various information, packaged similarly to the first image.  But I need just a single line (e.g., 0402), not all of them.

I'm ultimately trying to get the contour lines separate so I can create a DEM from this data.

Any suggestions?


texts: What's in the alpha channel in jp2 imagery file downloaded from usgs?

texts:
I'd like to know what's in the 4th channel in jp2, imagery files downloaded from USGS. In Photoshop those files open in grayscale mode, so to fix it I change them to multichannel and then RGB. But I assume it doesn't matter. From the beginning file contain 4 channels. RGB + something...

It looks like this:


texts: Sum integer values if values in another column are equal in QGIS?

texts:
I want a new column called "Summed" containing the summed value per "Populatie" name.

How can I accomplish this?


texts: Print Maps - OpenLayers, GeoServer, GeoExt, MapFish

texts:
Can anyone point me towards a clear and unambiguous description of how to install and code a very simple map print function for OpenLayers/GeoServer?

I am using GeoServer 2.5.2 and OpenLayers 2.12 on a Windows machine. Javascript handles the OpenLayers map and I send a lot of stuff to VB.NET for heavier backend work.

As of now I have added the GeoServer 2.5 printing plugin, instructions here:

http://docs.geoserver.org/2.5.x/en/user/community/printing/

So far so good, everything works. After this point there are two options, MapFish Print Module or GeoExt. It appears that GeoExt requires Extjs which I would have to license (my project is not opensource) and I don't want to deal with licensing or license fees. Can GeoExt be used alone with OpenLayers/GeoServer to print?

The other option is MapFish - MapFish Print Module - which seems like a much better alternative but I am completely befuddled on how to integrate it into my existing GeoServer/OpenLayers/JavScript/Vb.NET application. For instance it says to build the file. Build it where? In windows command line? In a java editor?


texts:
I want to change the width of line features.

Is it possible to set the parameter in map units (meter) and not in pixels?

I know that I can apply a buffer to get the same result. But I think it should be possible to do that directly.


texts: Exporting data from Oracle to Esri File Geodatabase grid size error

texts:
My specs: ArcInfo 10.1, Oracle 11g

I have a line dataset which is sat in Oracle 11g.  I can read/query/symbolise this in ArcDesktop without any problems.

I now want to export a selection into a file geodatabsae, but I keep getting this error:

The name it gives in the error (U1BLine) is the name I am trying to give my exported dataset.

I have tried exporting only the selection and the whole dataset.  I get the same error with both.

I have also tried to export it to a shapefile.  This works.

I know I could export to shapefile and then move that into a geodatabsae, but I'd like to be able to do it in a single step.

Other datasets from the same Oracle database have exported without a problem.

Both the dataset and the dataframe are set to British National Grid.

Anybody know how to fix this error?


texts: “addItem” in Portal for ArcGIS not working?

texts:
I'm trying to upload a PDF file to my content in Portal for ArcGIS using python. So far, I'm able to succesfully create a token, but I'm having troubles using REST operation :  Add Item/Add ItemMultipart/Add Part. Here's my code:

def get_token(username, password):
  """ Returns an authentication token for use in ArcGIS Online."""
  portalUrl = 'some_portal_url'
  token_url = 'some_portal_url/sharing/rest/generateToken?'

  # Set the username and password parameters before
  #  getting the token.
  params = {'username' : username,
            'password' : password,
            'client' : 'referer',
            'referer': portalUrl,
            'expiration': 1209600,
            'f' : 'json'}

  token_response = encode(params,token_url)

  if "token" in token_response:
      token = token_response.get("token")
      return token

def encode (params, url,):
  data = urllib.parse.urlencode(params)
  data_encoded = data.encode("utf-8")
  request = urllib.request.Request(url, data=data_encoded)
  context = ssl._create_unverified_context()
  response = urllib.request.urlopen(request,context= context)
  content = response.read()
  content_decoded = content.decode("utf-8")
  return json.loads(content_decoded)

def uploadFile(params, url):
  '''Creates a new item in a user's content.'''
  return encode(params, url)

sourceToken = get_token(username='username' ,password='password')
pdfFile = open('C:\\some\\pdf\\path.pdf')

params={'fileName': 'E','multipart' :'true','f' :'pjson','token':sourceToken}
addItemUrl='some_portal_url/sharing/rest/content/users/username/addItem?multipart=true&'
reply = uploadFile(params,addItemUrl)
itemID = reply["id"]

#From here it doesn't work

parameters = {'file' : pdfFile,  'partNum': '1','token' : sourceToken,'f' : 'pjson'}
url = 'some_portal_url/sharing/rest/content/users/username/items/'+itemID+'/addpart'
r1 = uploadFile(parameters,url)

parameters = {'token' : sourceToken,'f' : 'json'}
url = 'some_portal_url/sharing/rest/content/users/username/items/'+itemID+'/commit'
r2 = uploadFile(parameters,url)


texts: ArcGIS error when using hfactor table in Path Distance

texts:
I am trying to calculate isochrones for sailing under various weather conditions.

Inputs for Path Distance
Input raster or feature source data: .shp with coastal points
Input cost raster: Raster of wind speed
Input horizontal raster: Raster of wind direction (0-360 degrees in integer format)
Horizontal factor: ASCII Table of 0 to 180 with corresponding values

I can use any other horizontal factor (Binary, Forward, Linear, Inverse Linear, etc) and get a valid result, but when I use a table, I receive ERROR 999998: Unexpected Error. Failed to execute (PathDistance).

Anyone encountered this or found a workaround? I had to do a lot of troubleshooting, and at this point, I'm pretty confident it's something to do with using a table as the hfactor.


texts: Knowing the pixel size for a raster before it was resampled

texts:
Is there a method to know the exact pixel size (the original size) of a raster “aerial photo”?

I have purchased an Arial photo with 10 cm pixel size resolution, and I received a sample of it, but the issue here that I feel that the resolution on the image is not 10 cm , I come up with this conclusion by comparing these samples with other 10cm samples from other projects and the old aerial photos visually more clear than the new sample.

Both Erdas and ArcMAP shows the same pixel size (10cm) for both samples, and As I know you can produce any GSD in erads form the same Aerial photo by resampling the image to smaller pixels, and here the resolution stays the same but the pixels are smaller,

What I am trying to prove that this aerial photo was resampled form a higher GSD (maybe 15cm) to 10cm, is that possible?


texts: UserWarning when trying to get centroid from a polygon GeoPandas

texts:
I have this geopandas dataframe:
df1 = gpd.GeoDataFrame(df,crs=rasterio.crs.CRS.from_string("epsg:4326"), geometry=df['geojson'].apply(json_load))
df1 = df1.set_crs(epsg=4326)
print (df1.shape)
df1.head()

 account_id area    geojson                                             geometry
0   1365    1.7     {"type": "MultiPolygon","coordinates":[[[[35.4...   MULTIPOLYGON (((35.48632 32.78037, 35.48600 32...
1   3019    0.6     {"type": "MultiPolygon","coordinates":[[[[4.92...   MULTIPOLYGON (((4.92786 43.68305, 4.92951 43.6...

df1.crs
>> <Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

I would like to get the coordinates of the centroid of each polygon.
I used to be able to get it by df1['geometry'].centroid without any warning, but now when I try that line I get:
UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  """Entry point for launching an IPython kernel.

I don't understand why because my CRS looks fine, also the coordinates of the centroid look fine, should I just ignore this warning?
My packages versions:
print(pyproj.__version__)
print(rasterio.__version__)
print(gpd.__version__)

>> 2.6.1.post1
>> 1.1.5
>> 0.8.1


texts: Iteratively Change Object's Instance Attributes PyQGIS

texts:
I am trying to create new instance attributes (NOT feature attributes) for a collection of features in a vector layer in PyQGIS. I'm getting some weird results and am wondering if this is even possible. For instance, I can change the instance attributes of a single feature:

feat = iface.activeLayers().getFeatures().next()

feat.NEW = 'attribute'

And it seems to commit to that feature:

dir(feat)

['NEW', '__class__', '__delattr__', '__delitem__', '__dict__', '__doc__', '__format__', '__geo_interface__', '__getattribute__', '__getitem__', '__hash__', '__init__', '__iter__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setitem__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'attribute', 'attributes', 'deleteAttribute', 'fieldNameIndex', 'fields', 'geometry', 'geometryAndOwnership', 'id', 'initAttributes', 'isValid', 'setAttribute', 'setAttributes', 'setFeatureId', 'setFields', 'setGeometry', 'setGeometryAndOwnership', 'setValid']

feat.NEW

'attribute'

However, when I try to iterate through the features using a feature iterator, I get no such luck!

for feat in iface.activeLayer().getFeatures():
    feat.NEW = 'attribute'

Then check:

feat = iface.activeLayer().getFeatures().next()

feat.NEW

Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsFeature' object has no attribute 'NEW'

I'm assuming that the feature iterator creates some objects that have access to the features within the layer but does not actually point to the feature objects themselves... does anyone know of a way to iteratively change the instance attributes of PyQGIS objects? Is this even possible?


texts: GRASS GIS: how to define several variables in a loop without the command locals()?

texts:
I am working on a Python script in GRASS GIS and want define several variables in a loop. My simple loop is working, but I don't know how to integrate sevaral variables.

RG_koor=g.read_command("v.out.ascii", input="e_"+Nr_RG[i], fs=",", format="point")
zeil=RG_koor.split("\n")
easting=[];northing=[]  # define northing and wasting koordinates, empty

for i in range(len(zeil)-1):
  zeile=[]  # define "zeile" for easting and northing coodinates
  zeile=zeil[i].split(",") #split the colums when the comma appears
  easting.append(float(zeile[0])) #easting, all values before the first comma
  northing.append(float(zeile[1])) #northing, all values before the second comma

I tried the command locals (), but I read, that it is not recomended.Do you have  an idea, how I could improve my script without using the command locals()? I only would like to generate the northing and easting values of several data in two different files. The upcoming code is working until the second .split. So, I still do not understand, why the second .plit does not work.

Nr_RG=["001","002","003","004","005"]
for i in range (len(Nr_RG)):
    locals()["RG_koor{0}".format(i)] = g.read_command("v.out.ascii", input="e_"+Nr_RG[i], fs=",", format="point")
    locals()["zeil{0}".format(i)] = locals()["RG_koor{0}".format(i)].split("\n")
    locals()["easting{0}".format(i)]=[];locals()["northing{0}".format(i)]=[]  # define northing and wasting koordinates, empty

    for i in range (len(locals()["zeil{0}".format(i)])):
        locals()["zeile{0}".format(i)]=[] # define "zeile" for easting and northing coodinate
    locals()["zeile{0}".format(i)] = locals()["zeil{0}".format(i)].split(",") #split the colums, when the comma appears
    locals()["easting{0}".format(i)].append(float(locals()["zeile{0}".format(i)][0])) #easting, all values before the first comma
    locals()["northing{0}".format(i)].append(float(locals()["zeile{0}".format(i)][1])) #northing, all values before the seocnd comma

Here are my outputs for the second command:

    Nr_RG=["001","002","003","004","005"]
for i in range (len(Nr_RG)):
    sqlite="""substr(cmt,1,4)='R"""+str(Nr_RG[i]+"'")
    print sqlite
    g.run_command("v.extract", overwrite=True, input="rec_RG", output="e_"+Nr_RG[i], where=sqlite)
    locals()["RG_koor{0}".format(i)] = g.read_command("v.out.ascii", input="e_"+Nr_RG[i], fs=",", format="point")
    locals()["zeil{0}".format(i)] = locals()["RG_koor{0}".format(i)].split("\n")
    locals()["easting{0}".format(i)]=[];locals()["northing{0}".format(i)]=[]

substr(cmt,1,4)='R001'
0
substr(cmt,1,4)='R002'
0
substr(cmt,1,4)='R003'
0
substr(cmt,1,4)='R004'
0
substr(cmt,1,4)='R005'
0

When i split the second time by using locals (), the command is not working.

 for i in range (len(Nr_RG)):
    locals()["zeile{0}".format(i)]=[] # define "zeile" for easting and northing coodinate
    locals()["zeile{0}".format(i)] = locals()["zeil{0}".format(i)].split("\n")

Traceback (most recent call last):
  File "<input>", line 3, in <module>
AttributeError: 'list' object has no attribute 'split'


texts: How to load valid PostgreSQL QgsVectorLayer

texts:
Note:

This is NOT a closed question. The following links does NOT provide an
answer to the solution:
Add PostGIS layer to QGIS via Python Console
The Vector layer still remains invalid with the recommended changes
being made.
PyQGIS QgsVectorLayer() Loading Invalid Layer in Standalone Python Script?
I can also confirm my QgsApplication prefixpath is correct

I have a QGIS script I am writing in PyCharm, that I am trying to load a vector layer stored in a PostgreSQL database. When I print the layer's isValid() method I get False. Here is my code:
from qgis.core import *

db_client = 'postgres'
db_host = 'localhost'
db_port = '5432'
db_name = 'database'
db_user = 'user'
db_password = 'pass123'
db_schema = 'public'
tablename = 'Geo_Layer'
geometrycol = 'geom'
tract_number_index = 3

QgsApplication.setPrefixPath('/usr', True)
qgs = QgsApplication([], False)
# I have also tried the following:
# qgs = QgsApplication([], True)
qgs.initQgis()

geo_uri = QgsDataSourceUri()
geo_uri.setConnection(db_host, db_port, db_name, db_user, db_password)

geo_uri.setDataSource(db_schema, tablename, geometrycol, '', 'id')
# I have also tried the following
# geo_uri.setDataSource(db_schema, tablename, geometrycol)

geo_layer = QgsVectorLayer(geo_uri.uri(False), "Test", "postgres")
# Other configurations I have tried
# geo_layer = QgsVectorLayer(geo_uri.uri(), "Test", "postgres")
# geo_layer = QgsVectorLayer(geo_uri.uri(), "Test", "ogr")
# geo_layer = QgsVectorLayer(geo_uri.uri(False), "Test", "ogr")
# A working version using an exported geopackage
# geo_layer = QgsVectorLayer('/home/<path to geopackage>/geo_layer.gpkg|layername=geo_layer', "Test", "ogr")
print(geo_layer.isValid())

qgs.exitQgis()

I have provided the other QgsVectorLayer configurations I have tried. All print that the layer is not valid.
QGIS Version: 3.16.3-Hannover Python Version: 3.8.5 Ubuntu Version: 20.04.02 LTS
I have check my credentials with DBeaver and I am able to connect.
I have made sure I have the correct QgsApplication prefixPath by opening my QGIS Desktop, opening the Python Console and executing the following:
QgsApplication.prefixPath()
'/usr'

Also in QGIS Desktop the specific vector layer I am trying to validate, I can export that layer into a GeoPackage and connect to it and validate the Vector Layer. So I think I have the correct QgsApplication configuration.
geo_layer = QgsVectorLayer('/home/<path to geopackage>/geo_layer.gpkg|layername=geo_layer', "Test", "ogr")

This is the workaround I have in place, but I would like to connect directly to my PostgreSQL database.


texts: Getting AttributeError: 'NoneType' object has no attribute 'sourceCrs' from TIN Interpolation in QGIS?

texts:
I am trying to use the TIN interpolation to interpolate conductivities between a vector layer in QGIS 3.4. and I always get this error:
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\qgis\TinInterpolation.py", line 171, in processAlgorithm
crs = layer.sourceCrs()
AttributeError: 'NoneType' object has no attribute 'sourceCrs'

Does anyone have encountered this error already? I didn't find anything about it. The CRS is homogeneous in all my layers.


texts: Restoring or repairing non-working QGIS file (*.qgs)?

texts:
It is the second time in a week that have lost the information in the project file (.qgs). When I load the project the only informationI see is the groups in the layers panel, but no layers.

I have no idea what went wrong, inspecting the .qgs file I see that all my layers, coniguration layers, data sources and a lot more are there (20k lines !!). My guess is that something could have gone wrong with some virtual layers, so I think I could get ride of them, but how?

I wonder if there is a way to retrieve some data and info in order to setup the project again from the .qgs file to build a project brick by brick and not start from zero.

System information: Windows 7 Enterprise, QGIS version 2.18.14

Edited:

I tried to use the .qgs~, but the problem persisted


texts: Calculating length of curved line derived from intersection between two polygons using ArcPy?

texts:
I'm trying to get the density of intersection points along their related ellipse. What I'm missing is the length of the arc of the ellipses comprised in between the yellow polygon (which is a landslide deposit).
I would need help about step by step algorithm to calculate the distance between two points along a curved line.


texts: Deploy/ Copy Feature Class under users AD account with gis Schema,

texts:
This question is regarding to SQL Server 2017 and ArcMap or ArcMap Pro environment.
Little history: Currently GIS team deploying/Copying Feature class (FC) from ArcMap 10.4 or Pro under one SQL account and when they copy the FC under that account it by default creates Feature class with GIS Schema (e.g., gis.city, gis.metro etc) and in back end tables created with gis schema as well.
New Approach:
I am planning to use their AD account to deploy the feature class so that way one account password not to be shared and for more security and audit purpose as well.
When I am trying to Copy FC (from Old Srver\geodatabase to New Server\New Geodatabase) under AD account following things happens.
Even though my AD account has default Schema assigned gis not dbo, Still FC gets deployed under dbo Schema (e.g dbo.city) from ArcMap. And my account has db_owner permission.

When my AD account has Read/Write Permission and Db_Owner Removed I get below error. copying FC from old server to new server with different database name as well.
" Paste Failed. Failed to Paste gdb.gis.city database user name and current user schema do not match. [gdb1,"contoso\xxxx".city] DBMS table not found."
I am looking to for...

How multiple users can deploy/copy/import/ Edit FC under their AD account with keeping gis schema so by default it goes under gis.tablename in SQL always.
what database level permissions requires for AD users in sQL since I have been trying with various permissions and no success yet.
Is this a good practice then keeping one SQL account for 10-15 GIS users to perform their deploy/copy/import/Edit tasks?


texts: Joining attributes of census tracts with neighborhoods in ArcGIS when the boundaries don't match up perfectly?

texts:
I have a shapefile of Montreal census tract boundaries (below left), where for each tract I've joined data for 17 census variables including population, land area, total private dwellings, etc.

I also have a shapefile of Montreal neighborhood boundaries (above right) which are larger than the census tracts. My goal is to join the attributes of the layers so that I have a resulting layer of census tracts, where each tract now has a neighborhood name associated with it. I will then use this layer to aggregate data by neighborhood.

The big problem here is that the boundaries of the census tracts sort of align with the neighborhoods, but NOT perfectly. Using tools like intersect isn't working because the census tracts that are not aligned perfectly end up getting repeated in the final layer for every neighborhood that they cross into (because they don't neatly fit into a single neighborhood).

Is there anything I can do to solve these boundary issues or any other tool that could help with this?


texts: How to filter and show properties of features when clicked using Leaflet?

texts:
So, I'm mapping some features using Leaflet and want to be able to filter by date of the feature and also, when clicked, the feature will display some properties about it. I have yet to implement the date portion so I'm using a static date to filter through my data set and I'm struggling with getting each point to display information about it when clicked.
All other solutions that I've looked through use the portion where I've implemented the filter but I need that.

        <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css"
        integrity="sha512-xodZBNTC5n17Xt2atTPuE1HxjVMSvLVW9ocqUKLsCC5CXdbqCmblAshOMAS6/keqq/sMZMZ19scR4PsZChSR7A=="
        crossorigin=""/>

I've gotten all of the points to display "Hello" when clicked using:
filtered_points.bindPopup("Hello);

But I haven't been able to get them to display properties about each respectively.


texts:
i'm trying to layout a map in ArcGIS09, showing samplings in a specific region:
i have a layer with samplings on different sampling dates and sites. sometimes, there are several samples at one site, done at different dates.

the attribute i want to layout is the population structure (very good, good, moderate, poor, bad), which is easy.
but now, i would like to sort the points in a way, that the older sampling dates are hidden by the more recent ones. and it would be even nicer, to layout the older sampling dates bigger than the more recent one, in a way that they are still visible below...

is this possible in any way?

thanks a lot for your help!


texts: Is there anyway you can invert the fill of a single polygon on QGIS when you have several polygons in the layer?

texts: I have used QGIS inverted polygon fill too but was wondering if you could isolate individual polygons - a bit like the atlas feature?


texts: No "committedFeaturesAdded" signal after dataProvider().addFeatures() - QGiS 2.14.8

texts:
I'm working on plugin updating existing shapefile by adding new features to it. And I have a problem with signal that hasn't been emitted by QGIS, when I've added features to layer, using dataProvider().addFeatures().
I'm pretty sure that I have defined correctly connection signal and function, when it is emitted.

layer.committedFeaturesAdded.connect(self.onFeaturesAdded)

Becasue when I add feature using tool in QGIS from Digitizing Toolbar, it works fine and function is running.

The problem is when I add features manually - using python.

def selectFeatures(self,diviLayer, userLayer84):
        self.diviLayer=diviLayer
        self.userLayer84=userLayer84
        features=[]

        Processing.initialize()

        general.runalg('qgis:selectbylocation', self.userLayer84, self.diviLayer, u'intersects',0.000001,0)
        self.userLayer84.invertSelection()
        number=self.userLayer84.selectedFeatureCount()

        for feature in self.userLayer84.selectedFeatures():
            attributes=[feature['Lokaliz'],feature['Zadanie'],feature['Projekt'],feature['Wykonawca'],feature['X'],feature['Y'],feature['h'],feature['Azymut'],feature['Dokl'],feature['KmEt'],feature['Km2Et'],feature['Poloz'],feature['Domiar'],feature['Nr_dz'],feature['Obr'],feature['JE'],feature['Pow'],feature['Woj'],feature['Data'],feature['Czas'],feature['xtrPhtC'],feature['xtrComC'],feature['xtrAudC'],feature['xtrFilC'],feature['NazwaRap'],feature['Sciezka']]
            feat=QgsFeature()
            feat.setAttributes(attributes)
            feat.setGeometry(QgsGeometry.fromPoint(QgsPoint(float(feature['Y']),float(feature['X']))))
            features.append(feat)

        return [features,number]

    [features,number]=self.selectFeatures(self.diviLayer,self.userLayer84)

    self.diviLayer.startEditing()
    self.diviLayer.dataProvider().addFeatures(features)
    self.diviLayer.commitChanges()

It seems like the signal is not emitted? Has anyone of you had a similar problem? Or maybe you know, where I'm making mistake.


texts:
I have some substances TPH, CO2, SO2, VOC, analyzed in a period of time 2001 -  2017. In the .mxd each year is a group layer, containing all the substances for that year:

2001:
      - TPH

  - CO2

  - SO2

  - VOC

2002
      - TPH

  - CO2

  - SO2

  - VOC

...until 2017

I want to apply the same symbology for TPH substance in 2001, 2002, 2003 etc. I have made .lyr files for all the substances, but I have managed to write a functional code to apply one symbology for all layers in the same group layer.

Any solution to apply the same symbology for just one layer in all the group layers?

I have many more substances than the ones mentioned!

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")

df = mxd.activeDataFrame

groupName = "2001"

allLayers = arcpy.mapping.ListLayers(mxd, groupName, df) [0]

for layer in allLayers:

    arcpy.ApplySymbologyFromLayer_management(layer, "TPH")

PS: I use ArcMap 10.6 and the batch function for the Apply Symbology From Layer tool is not working


texts:
I'm new to Python and trying to work with TIFF images using the following  code:

from PIL import Image
import numpy

im = Image.open('LakeEx.tif')
imarray = numpy.array(im)
imarray.shape
im.show()

and I got this error:

  "ValueError: tile cannot extend outside image"

Is there any specific criteria to the TIFF file to work with Python: size, bit,... etc?


texts: Add Geodjango MultiPolygon to Leaflet Map

texts:
I'm building a geospatial mapping application with Geodjango. I have created a model which has a MultiPolygonField():
class newJob(models.Model):
    job_name = models.CharField(max_length=64)
    job_user = models.CharField(max_length=64)
    job_desc = models.CharField(max_length=64)
    job_loc = models.MultiPolygonField()

I have a simple function-based view to gather the required data from the database:
def viewjob(request):

        req = request.GET
        print(req)

        job_name = newJob.objects.values_list('job_name', flat=True).get(pk=req['search'])
        job_desc = newJob.objects.values_list('job_desc', flat=True).get(pk=req['search'])
        job_loc = newJob.objects.values_list('job_loc', flat=True).get(pk=req['search'])

        print(job_loc)
        print(type(job_loc))

        context = {

        'job_name': job_name,
        'job_desc': job_desc,
        'job_loc': job_loc

        }

        return render(request,'ecomap/viewjob.html', context)

The variable job_loc is of class django.contrib.gis.geos.collections.MultiPolygon. An example looks like:
SRID=4326;MULTIPOLYGON (((0.9008789062500006 52.42252295423906, 0.3900146484375003 52.05586831074773, 0.5657958984375002 51.73043054294018, 1.2469482421875 51.90700188674099, 1.6754150390625 52.22443459871997, 1.2139892578125 52.43926935464697, 0.9008789062500006 52.42252295423906)))

My basic template looks like:
{% extends "base.html" %}

{% block body %}

{% load leaflet_tags %}
{% leaflet_js %}
{% leaflet_css %}

{{ job_name }}

{{ job_desc }}

  {{job_loc}}

Your Area

  {% leaflet_map "yourmap" callback="window.map_init_basic" %}

 {% endblock %}

I was hoping it would be as easy as passing job_loc to Leaflet given the data is spatial but I'm thinking that there are perhaps more intermediate stages to get the data in the right format for plotting.
What is the easiest way to plot the polygon represented in job_loc on a leaflet map embedded in my template?


texts: GDAL GeoTIFF to GeoPackage overview building failed error

texts:
I have multiple GeoTIFFs in UTM projection. I need to create GPKG SQLite databases. I noticed that some (Example3) conversions are failing with overview error. What is cause for the failure ? How could we overcome it?
Example1:
D:\Raster_Data\temp>gdalwarp -t_srs EPSG:4326  Kan_UTM.tif Kan_WGS84.tif
Copying color table from Kan_UTM.tif to new file.
Creating output file that is 18361P x 10657L.
Processing Kan_UTM.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GTiff -a_srs EPSG:4326 -expand rgba Kan_WGS84.tif Kan_WGS84_RGBA.tif
Input file size is 18361, 10657
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GPKG -ot Byte -co TILING_SCHEME=GoogleMapsCompatible Kan_WGS84_RGBA.tif Kan.gpkg -a_srs EPSG:4326
Input file size is 18361, 10657
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdaladdo Kan.gpkg 2 4 8 16 32 64 128 256 512 1024 2048
0...10...20...30...40...50...60...70...80...90...100 - done.

Example2:
D:\Raster_Data\temp>
D:\Raster_Data\temp>gdalwarp -t_srs EPSG:4326  Pho_UTM.tif Pho_WGS84.tif
Copying color table from Pho_UTM.tif to new file.
Creating output file that is 19376P x 11003L.
Processing Pho_UTM.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GTiff -a_srs EPSG:4326 -expand rgba Pho_WGS84.tif Pho_WGS84_RGBA.tif
Input file size is 19376, 11003
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GPKG -ot Byte -co TILING_SCHEME=GoogleMapsCompatible Pho_WGS84_RGBA.tif Pho.gpkg -a_srs EPSG:4326
Input file size is 19376, 11003
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdaladdo Pho.gpkg 2 4 8 16 32 64 128 256 512 1024 2048
0...10...20...30...40...50...60...70...80...90...100 - done.

Example3:
D:\Raster_Data\temp>gdalwarp -t_srs EPSG:4326  Chic_UTM.tif Chic_WGS84.tif
Copying color table from Chic_UTM.tif to new file.
Creating output file that is 20121P x 10164L.
Processing Chic_UTM.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GTiff -a_srs EPSG:4326 -expand rgba Chic_WGS84.tif Chic_WGS84_RGBA.tif
Input file size is 20121, 10164
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdal_translate -of GPKG -ot Byte -co TILING_SCHEME=GoogleMapsCompatible Chic_WGS84_RGBA.tif Chic.gpkg -a_srs EPSG:4326
Input file size is 20121, 10164
0...10...20...30...40...50...60...70...80...90...100 - done.

D:\Raster_Data\temp>gdaladdo Chic.gpkg 2 4 8 16 32 64 128 256 512 1024 2048
Warning 1: Use of overview factor 512 causes gpkg_zoom_other extension to be needed
ERROR 1: Could not find dataset corresponding to ov factor 512
Overview building failed.


texts: Why is the reciprocal heading of a line drawn between two points on Google Earth not equal to 180 degrees?

texts: If I measure a line from point (A) to point (B) Google Earth gives me a heading. if I measure a line from B to A I would think it should be 180 degrees from (A) to (B). It is not. I'm trying to find the angles of a triangle drawn on Google Earth and the variation is making it difficult.


texts: Zooming to feature selected from QListWidget using PyQGIS

texts:
I have displayed in a QListWidget several features ranked according to their values from a field named "CODE". Based on this list I would like to select a feature and after a click on a QpushButton to zoom to the extent of this feature in the canvas.
Here is my crude attempt (I am still a beginner):
def append_features_to_QlistWidget (self): #associated to QpushButton_1

    #----------------function sample--------------------

    self.layer = QgsProject.instance().mapLayersByName('my_layer')[0]

    self.selection = []

    for self.feat in self.layer.getFeatures():
        self.attrs = self.feat.attributes()
        self.code= self.attrs[0]
        self.selection.append(self.code)

    self.dlg.listWidget.addItems(self.selection)

def zoom (self): #associated to QpushButton_2

    inds = [i.text() for i in self.dlg.listWidget.selectedItems()]

    attribute_value = {}

    for i in inds:
        attribute_value[i] = [f for f in self.layer.getFeatures() if f["code"] == i][0]

    self.layer.selectByIds(attribute_value)

    iface.mapCanvas().zoomToSelected(self.layer)

    ##Other incomplete attempt:
    # ids = [i.id() for i in self.selection]
    # self.layer.selectByIds(ids)
    #iface.mapCanvas().zoomToSelected(self.layer)


texts:
I'm trying to download tiled Daymet climate data ([see here])1. One can automate this process using a shell script, but you need to know the tileIDs that you want.
Is there a way to download the tile grid as a vector file, so I can easily get a list of tileIDs using a spatial selection? Or any other way of selecting tileIDs within geographic bound?


texts: Linking SQL Server database to ArcGIS Online

texts:
I'm not sure this is possible, and the ESRI documentation is a bit of a mess it seems. Or maybe it's not possible at all and the ESRI documentation is giving me the nearest it can.

Basically, the organisation I work for has a fair bit of data on an SQL Server - both spatial and non-spatial.

We want to link this to ArcGIS Online / Collector, so that a user can log in and see the latest data and have it updated. Bonus points if they can make updates back.

I believe this will require ArcGIS Server, and there seem to be a whole myriad of options. I've been told that it doesn't, but not sure whether to believe it. Also, we'll be trialing ArcGIS Server shortly, in order to make our SQL Server database into a geodatabase (this is the approach ESRI themselves have said to use).

If we do that, even once the trial expires, would we be able to have ArcGIS Online look at that?

As near as I can understand it, we need to set up an ArcGIS Server web service.

Would we still be able to use this once the trial expires?

We've currently got a standard editor licence.


texts: I receive the following non descriptive "error reported on layer ...." when trying to checkout a feature class from a PGDB for disconnected editing.  All other feature classes checkout as background shapefiles from the same database but will get the same error if I try them as an AXF file.  ArcPad 10.2.3, ArcGIS 10.1 and Windows 10.


texts: QGIS says "osgeo [python-gdal]" isn't installed on Mac?

texts:
I know similar questions have been asked (like Missing QGIS raster menu tools - '"osgeo [python-gdal]" module is missing') -- but I'm encountering a unique problem. I installed GDAL complete 1.11 from KyngChaos, but when I open QGIS, I get the following error message:

Interestingly enough, when I look at the Libraries directory in Finder, there isn't a subfolder called Frameworks. Which tells me something went wrong.

I'm running Mac OS X 10.9.2.


texts: How to join (draw new) lines that are close but not touching - large dataset

texts:
As per the picture below I've got 2 lines that are very close but not touching (0.013mm) in a T-junction format but could be any kind of intersection.

      |                           |                   |
----  |   needs to become   ------|    ...or    ------|--
      |                           |                   |

I don't know how many there are of these (10k features in total in this dataset) but it's breaking a path finding algorithm so need them to magically get joined together.

A new line segment is fine but moving the existing one would be better (as as not to increase the existing number of features unnecessarily).

I looked at snapping but it appears that would then be a manual operation which is not ideal (as there are potentially multiple datasets all with 10K plus features) though could work if I can generate a list of these 'near but not touching' lines and there are not many of them.

I'm wondering if there's any buffering function that could do it but then not sure how to move the existing node or create the new line based on that.

Any thoughts on how to do this as a bulk operation in QGIS (3.4.11)

Thanks in advance!!


texts: How to Sample Locations from a Satellite Image

texts:
I am using the Near infrared and red bands from satellite images to research soil quality in southern Africa. The study site is about 20 km wide and 20 km long so almost square shape, and I am using 5 m resolution satellite imagery.

What is the best way to sample different locations so that I get enough samples and also to minimize spatial autocorrelation?

Is grid sampling method the best?

I am new to studying soil.


texts: When drawing polyline using QGIS it starts on outside of cursor and not in centre?

texts:
When drawing a polyline in QGIS 3 it starts on the outside of the cursor, not in the centre. Same thing when measuring a line. This makes accuracy difficult.

Am I missing something?

I am using OS X (High Sierra).


texts: Importing data from Google Earth (points, lines and polygons with their attributes) into QGIS?

texts:
I want to import data in QGIS from Google Earth (points, lines and polygons) similar to OSM.

Is it possible to do that and, if yes, then how?


texts: PostGIS runs out of memory due to the use of spatial index

texts:
I am running this query on PostgreSQL.

CREATE TABLE table3 AS (SELECT
    t1.id
FROM
    table1 AS t1
    CROSS JOIN LATERAL
        (SELECT
            t2.geometry AS geom_closest
        FROM table2 t2
        ORDER BY t1.geometry <-> t2.geometry
        LIMIT 1) AS nn);

To explain it quickly: I have a grid of points in table1, a list of features in table2. I want to find, for each point in table1, the minimum closest feature in table2 (for example, say I have in table2 the roads in the state of NY, and in table1 a regular grid over the state of NY, then for each point in my grid I want to find the closest road).

I have 8GB of RAM on my server. My tables are pretty small all things considered.

 tablename                  |    indexname                   |  num_rows   | table_size | index_size | unique | number_of_scans | tuples_read | tuples_fetched
table1                      | table1_idx                     | 1.47763e+06 | 96 MB      | 82 MB      | N      |               0 |           0 |              0
table1                      | table1_pkey                    | 1.47763e+06 | 96 MB      | 32 MB      | Y      |               0 |           0 |              0
table2                      | table2_pkey                    |        2916 | 8640 kB    | 80 kB      | Y      |               1 |        2916 |              0

Now, if I create a spatial index on table2 (GIST index), my system go crazy with my query. The RAM used takes off and eventually (< 5 minutes) my server runs out of memory and crashes hard:

WARNING:  terminating connection because of crash of another server process
DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, because another server process exited abnormally and possibly corrupted shared memory.
HINT:  In a moment you should be able to reconnect to the database and repeat your command.
server closed the connection unexpectedly
    This probably means the server terminated abnormally
    before or while processing the request.
The connection to the server was lost. Attempting reset: Failed.

If I do not use a spatial index, then my memory usage stays at 0.5% and my query just takes forever.

What is wrong with my query? How can I improve it? I will need to run the same kind of query on much larger datasets and I do not understand why using my spatial index breaks everything.

postgis_version: 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1


texts: Displaying vertex coordinates of a polygon or line without creating a new layer

texts:
I am wondering if there is a way to display coordinates for each vertex of a polygon by only using expressions, without creating a new layer. So this question is basically inspired by Displaying vertex coordinates of a polygon on a map in QGIS and enhancing it.
I thought about creating an array containing all vertices as geometry. So I have created this expression:
array_foreach( -- Go throuth every element of the following array
    string_to_array(geom_to_wkt(nodes_to_points($geometry)),','), -- Create the array of WKT information splitted by comma [ 'MultiPoint ((0.2974026 51.8987013)', '(0.00909091 51.50… ]
        make_point( -- Create a point geometry for every array content
            to_real(if(
                length(regexp_substr( -- Only create X coordinate if
                    regexp_substr(@element,'\\d.*\\d.[0-9]'),'([^ ]+)'))>0, -- everything between the Brackets of 'MultiPoint ((0.2974026 51.8987013)' before the whitespace '0.2974026' is not empty like ''
                    to_real(regexp_substr(regexp_substr(@element,'\\d.*\\d.[0-9]'),'([^ ]+)')), -- Extract the X coordinate of 'MultiPoint ((0.2974026 51.8987013)' --> 0.2974026
                    NULL -- If invalid coordinate use NULL, maybe replace with 0
            )),
            to_real(if(
                length(regexp_substr( -- Only create Y coordinate if
                    regexp_substr(@element,'\\d.*\\d.[0-9]'),'(?<=\\s).*'))>0, --  everything between the Brackets of 'MultiPoint ((0.2974026 51.8987013)' after the whitespace '51.8987013' is not empty like ''
                    to_real(regexp_substr(regexp_substr(@element,'\\d.*\\d.[0-9]'),'(?<=\\s).*')), -- Extract the Y coordinate of 'MultiPoint ((0.2974026 51.8987013)' --> 51.8987013
                    NULL -- If invalid coordinate use NULL, maybe replace with 0
            ))
        )
)

Which returns an array like [ <geometry: Point>, <geometry: Point>, <geometry: Point>, <geometry: Point>, ... ], which represents all the vertex of the polygon. One can easily extract X and Y coordinates out of it by adding x() or y() in front of make_point(), or by just adding another array_foreach() like array_foreach(*expression above*,x(@element)).
How could I use such an array to label every vertex with its coordinates? I thought of extracting them from the array, but cant iterate over it. Would this be possible? Or something similar? I am kind of stuck now...


texts:
I am making an addin to create subsets and export them to a new workspace. I have stepped through my program and getting error at this particular line while running the addin:

pExportOp.ExportFeatureClass((IDatasetName)pFeatureClassName,
        pQueryFilter, pSelSet, pGeometryDef, (IFeatureClassName)pOutDatasetName, 0);

Can someone help that what the error could be.


texts: Most of the tutorial videos, they simply add those training samples. Looking at the attribute table they seems like, they in some specific format. Like all the training samples belonging to perticular class have one object id. However, if you create training samples using ArcGIS each polygon will have different object ID.


texts: Converting pre-ADF coverages to current shapefiles?

texts:
I have a set of old coverage files, published at 1995-1996 in a CDROM by an official mapping authority. They were, therefore, made with ArcView 1.x and don't have the ".adf" extension. The files that come into the single folder are: ARC, ARX, CNT, CNX, DBLBND, DBLTIC, LAB, LOG, PAL, PAR, PAX, TXT and TXX (all without any exntension), besides a PAT.DBF.

I've tried to import that into ArcGIS 10, without success. I've tried to open it in ArcView 3.3 by adding it as a New Theme, but nothing is shown as well. I've tried to convert them with ogr2ogr and with ArcINFO2Shape utility (deliverd at: http://freegeographytools.com/2007/converting-e00-vector-data-to-shapefiles-a-free-and-fairly-painless-approach). Nothing seems to work. I also don't have access to a an older ArcWorstation.

Is there any way to do achieve such a task?

PS: I am aware of the other two questions asking similar things here at the GIS-Stackexchange. My question, however, deals with even older non-adf extension data and, as described, the processes being discussed in the other related questions did not apply to my data.


texts: Geocoding USA addresses that cannot be sent over internet?

texts: For a one-time project, I need to geocode a few thousand addresses.  In the past I've used various online resources for this kind of thing (e.g., Google Maps API), but the addresses I'm working with have to be kept confidential - which means no sending it over the Internet, unless there's some iron-clad guarantee of privacy.  What other options do I have?


texts:
Is it possible to split up the individual values of a raster to make each value become a different layer?

I am aiming to create an animated map with TimeManager. However, the data I need for each layer of the animated map is within one raster.

The raster is split into 16 values from 2000-2016. Is there a way to take out each value to create a new layer?

Here is a picture of the raster for more clarity.


texts: Calculating distance between points crossing antimeridian using R?

texts:
How can I calculate the distance between two points (lon/lat) if any of the points exceeds the antimeridian (+-180 longitude)?

I am using the distVincentyEllipsoid {geosphere} to calculate the shortest distance between two points. But if any of the points crosses the antimeridian the distance will be very close to zero.

I have seen Intersection algorithm that correctly handles the 180° meridian and poles? but it does not use R.


texts: Seperate each row of points from shapefile with latitude and longitude information

texts:
I have a shapefile which has point geometry points with latitude and longitude values.The image when opened in QGIS looks similar to the below picture and I want to know whether I can separate each row of points using latitude and longitude information using GeoPandas or any other library in Python. I am new to GeoPandas and geo-processing.

My end goal is to find if there are any missing points in each row. The above image doesn't have missing points in each row but the data which I have has missing points in some of the rows.

Things that I tried so far:
I started sorting the points based on latitude and longitude information for finding a pattern to separate each row but it doesn't help me that much, so I would like to know if there are any other ways I can separate each row of data points.


texts: Implement ArcGIS Custom Features with C#?

texts:
Is it possible to implement custom features in .NET?

If so, how?

I suspect the issue is being able to do COM aggregation with .NET, but I can't find anything on that ...

Does anyone know if it is possible to create a runtime COM callable wrapper for a .NET class that supports COM_INTERFACE_ENTRY_AGGREGATE_BLIND ?


texts: Proxy issues using WFS from Cadcorp GeognoSIS

texts:
I have web and WMS/WFS servers hosted on the same machine, although on different ports.
I'm currently using OpenLayers to consume the WMS feeds which works well, but now I need to allow the user to extract info about objects so am trying to incorporate WFS into the OpenLayers setup.
Being on different ports, I've set up proxy.asp to do the WFS requests:

<%
set objHttp = Server.CreateObject("Msxml2.ServerXMLHTTP")
strURL = Request("url")
objHttp.open "GET", strURL, False
objHttp.Send

If objHttp.status = 200 Then
    Response.Expires = 90
    Response.ContentType = "application/xml;charset=UTF-8"
    Response.BinaryWrite objHttp.responseBody
    set objHttp = Nothing
End If
%>

I've tested the proxy using this URL:
http:///maps/proxy.asp?url=http%3A//:4326/OgcService/WFS.impl?service=wfs?request=getCapabilities

This returns the GetCapabilities which is identical to that returned if I access it directly, but without the proxy.

However, when using this WFS, via proxy, I can't get OpenLayers to display any content. At this stage it could well be my OpenLayers code, which I can post if necessary, but the same WFS URL via proxy doesn't work in MapInfo or Cadcorp MapModeller (it does with just the base URL). I suspect the proxy is the issue.

Any ideas what the issue is likely to be?


texts: Updating Field Based on another Field using UpdateCursor?

texts:
I have a feature class in a FGDB called LeasePlans. It contains a field called Property_Ref. I then have a CSV file with a field called SAPNO. For the feature class, I need to check does the value in Property_Ref occur in the SAPNO of the CSV and if so I need to add the word Archived to the new field called Status. Below is the code I have thus far, I'm not sure what is wrong with it, I'm getting an error ValueError: too many values to unpack but I suspect that is the least of my problems.

#Archives Historic SAP Numbers
import arcpy

fc = "U:/Data/Leases.gdb/LeasePlans"
changes = "U:/Data/To_be_Deleted.csv"

#New Field Criteria
field_Name = "Status"
field_Nametype = "Text"
field_NameVal = arcpy.ValidateFieldName(field_Name)

#Adds the new field if it doesn't already exist
fList = arcpy.ListFields(fc,field_Name)
if not fList:
    arcpy.AddField_management(fc, field_Name, field_Nametype, "", "", "")

#List comprehension for IDs
ids = [i for i, in arcpy.da.SearchCursor
       (changes, "SAPNO") if i != None]

Relevant_Fields = ['Property_Ref', 'Status']

#update cursor to Insert 'Archived' rows
with arcpy.da.UpdateCursor (fc, Relevant_Fields) as cursor:
    for SAPNO, in cursor:
        if SAPNO in ids:
            row[1] = 'Archived'
            cursor.UpdateRow()

del cursor


texts: How to export training area information?

texts: I did a supervised classification of a whole Landsat 8 image with 700 Rois (training areas) with the Semi Automatic Classification Plugin. I want to export the class names, the coordinates and the pixel size from the training areas.Furthermore I need spectral information from each pixel from all bands 2-7. I want to produce a excel file with all mentioned values.  The final data shall be used with SVM (Support Vector Machine) to improve my classification. I don't know how to export all of the values.


texts: Connecting two points in ArcGIS for Desktop?

texts:
I have two types of point shapefiles. These two points are linked with another attribute table. Each Point from column A is linked to single/multiple values in column B. Each of these object_IDs have their respective coordinates in another table.
eg.

 object_ID1   object_ID2
      1          3
      1          5
      1          9
      2          1
      3          7

where, each of these values represents the object id of the two different shape files. How to connect these two points with

a line (eg. plot a line between 1 and 3, 1 and 5, 1 and 9)
select the other value (eg. if I select one and run a python script (say), automatically select the other point/s)

The main objective is to visually link these two points anyway.


texts:
Using WMS I do a GetCapabilities request on a GeoServer instance.
This works in FireFox, but not in Chrome, Safari, and Internet Explorer.

The point is that in the web.xml from the server I did add:

"Access-Control-Allow-Origin: *"

So the XML output does not have anything which will put the Access-Control-Allow-Origin header in.

And the case is that in any other request this not an issue, especially when I use jsonp in the jquery AJAX call.

However I didn't find any way to get the Capabilities in a JSON object.
In other words I am stuck with XML and this doesn't seem to work.

Is there a way to do the request on jsonp?
If not how do I set-up the server so that the GetCapabilities XML does use the Access-Control-Allow-Origin?


texts: How to display 1,000,000 fishnet polygons in ArcMap?

texts: I am working on ArcMap. I have created a fishnet with 1,000,000 polygons for pixels of size 0.125 degree. I would like to display the map for a variable that varies across cells. I tried symbology and also converting fishnet to a raster data, but both takes forever. Is there any easier way to do it?