texts:
I have been writing a python script to convert an ascii grid into a geotiff file and when running the Create method in gdal i am receiving the following error message:

william@will-HP-EliteBook-850-G1:~/Smells$ python SmellsGeoTag.py ./SmellsTestData/test3/input/tile_1 ./output
Traceback (most recent call last):
  File "SmellsGeoTag.py", line 47, in <module>
    dst_ds = driver.Create(out , ncols, nrows, 1, gdal.GDT_Float32 )
  File "/usr/local/lib/python2.7/dist-packages/GDAL-1.11.0-py2.7-linux-x86_64.egg/osgeo/gdal.py", line 394, in Create
    return _gdal.Driver_Create(self, *args, **kwargs)
TypeError: in method 'Driver_Create', argument 3 of type 'int'

however unlike the problem here my ncols and nrows values are non zero ints, so i am rather stumped as to what the problem is

source code:

#!/usr/bin/python
import sys, getopt, os, numpy
from osgeo import osr, gdal

if len(sys.argv) != 3:
    print ("Please use the following arguements: input path, output path");
    exit();

inRoot = sys.argv[1]
outRoot = sys.argv[2]

#finding valid files
potentialDates = os.listdir(inRoot)
dates = []
for i in potentialDates:
    if os.path.isdir(os.path.join(inRoot, i)):
        dates.append(i)

for i in dates:
    potentialFiles = []
    files = []
    potentialFiles = os.listdir(os.path.join(inRoot, i))
    for x in potentialFiles:
        if x.endswith(".asc"):
            files.append(x)
    #processing each valid file found
    for y in files:
        with open(os.path.join(inRoot, i, y), "r") as f:
            linelist = [line for line in f]
            junk, xllcorner = linelist.pop(0).split()
            float(xllcorner)
            junk, yllcorner = linelist.pop(0).split()
            float(yllcorner)
            junk, ncols = linelist.pop(0).split()
            int(ncols)
            junk, nrows = linelist.pop(0).split()
            int(nrows)
            junk, cellsize = linelist.pop(0).split()
            float(cellsize)
            junk, nodataval = linelist.pop(0).split()
            float(nodataval)
            data = [[float(digit) for digit in line.split()]for line in linelist]
            driver = gdal.GetDriverByName("GTiff")
            if not os.path.exists(os.path.join(outRoot,i)):
                os.makedirs(os.path.join(outRoot,i))
            out = os.path.join(outRoot,i,y+".tif")
            dst_ds = driver.Create(out , ncols, nrows, 1, gdal.GDT_Float32 )
            dst_ds.SetGeoTransform( [xllcorner,cellsize,0,yllcorner,0,cellsize] )
            # set the reference info
            srs = osr.SpatialReference()
            srs.SetWellKnownGeogCS("WGS84")
            dst_ds.SetProjection( srs.ExportToWkt() )
            # write the band
            dst_ds.GetRasterBand(1).WriteArray(data)

exit()


texts: ArcGIS ClassBreaksRenderer isMaxInclusive doesn't affect JSON

texts:
In using ArcGIS's ClassBreaksRenderer, I've made use of setMaxInclusive(false) to allow the breaks to not include the maximum value. This works fine when using the renderer object itself (such as in a Feature Layer). However, when I convert the ClassBreaksRenderer to JSON using the toJson function, it does not change the resulting JSON object. I'm using the renderer in a dynamic layer (setting it through LayerDrawingOptions), and when I do that it is being treated as if that property were not set, which affects the resulting image. Looking at the network traffic, I see that it's sending the JSON for the renderer. I've included sample code below based on the ArcGIS renderer_class_breaks sample to give an idea of what I'm talking about.

My question is: is there a way to get JSON for a renderer that includes a reference to the isMaxInclusive property, or that behaves in the same way? Or is this just something that was not implemented in JSON?

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">

    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title>Class Breaks Renderer</title>

    <link rel="stylesheet" href="https://js.arcgis.com/3.19/esri/css/esri.css">


texts: Joining attributes by location for more than one located feature using QGIS

texts:
In the standard join attributes by locationfeature in QGIS 2.0/2.18, there are two options:

Take attributes first located feature
Take summary of intersecting features

Is there a function in QGIS that allows to take attributes of more than located feature?
That is, to get ALL qualitative values the polygons that intersect/overlap. With an attribute column for the value of the first located feature, another for the second located feature, etc. This should work in cases where located features can overlap (e.g. species habitat polygons where more than one species can occupy the same location) but the original polygon onto which we want attributes of the second layer has spatially separate polygons (e.g. watersheds from which we want to know which species are present).
In my case here: I am working with one polygon shapefile from which I want to obtain the attributes for another polygon shapefile. According to the output from the summary of intersecting features, there are a total of up to seven intersecting features for a given polygon in the file I want to join to.
It seems that it is not possible in QGIS given the comments but it would be good to confirm. I may make a feature request.


texts: Synchronizing on-prem PostgreSQL raster database with AWS RDS

texts:
I have a folder of rasters that updates daily that I load into a PostgreSQL database using the raster2pgsql command line tool. I want to make this database outward facing so that I can connect to it. I like the idea of storing this database on AWS, but I don't know if it's possible to

Add rasters directly to the AWS PostgreSQL database, or
Sync my on-prem database with the database on AWS.

I really would like to avoid having to re-import the entire database daily as I only need to keep adding one raster per day. Is this possible (and how would I go about doing this?) or is there a better tool for doing this?


texts: Combining between CON and Focal function in ArcGIS raster calculator?

texts:
How do you write this code in the raster calculator for ArcGIS 10.2?

CON(FOCALMIN([Burning], WEDGE, 4, 285, 345 ) == -1 , -1, [Burning])


texts: Creating clipped DEM and NDVI rasters with exact same numbers of rows and columns?

texts:
I have some large NDVI and DEM images of Kerala, India. I also have a shapefile of the subdistricts in Kerala.  I need to clip both the NDVI and the DEM so that there is an individual raster of NDVI and DEM for each of 14 subdistricts.

I used the clip function of the image analysis window to do this quickly, but the NDVI and DEM for each subdistrict don't have exactly the same numbers of rows and columns, which is a requirement for my application ( they will be converted to ASCII ).  They are usually off by one column or one row.

I don't understand how this can happen because I am using the exact same shapefile to clip both rasters.  This is in ArcMap.  Sorry if this isn't enough info or a newbie question, I'm still learning.

The cell size of the two rasters must also be the same, so I can't just export the rasters and force a different number of columns because that would change the cell size.


texts: Polygonizing OpenStreetMap MultiPolygon with Shapely?

texts:
I'm trying to import a very basic OpenStreetMap multipolygon made out of two outer lines into Shapely:

import json
from shapely.ops import transform
from shapely.geometry import mapping,shape
from shapely.geometry.polygon import Polygon

way = """
{
   "features":[
      {
         "geometry":{
            "coordinates":[
               [
                  [
                     [
                        143.1305688,
                        72.1178189
                     ],
                     [
                        142.884063,
                        72.1186622
                     ],
                     [
                        142.8833764,
                        72.1820139
                     ],
                     [
                        142.8833764,
                        72.1820139
                     ],
                     [
                        143.122329,
                        72.1818038
                     ],
                     [
                        143.1305688,
                        72.1178189
                     ]
                  ]
               ]
            ],
            "type":"MultiPolygon"
         },
         "properties":{
            "type":"multipolygon"
         },
         "type":"Feature"
      }
   ],
   "type":"FeatureCollection"
}
"""

json = json.loads(way)

print(json['features'][0]['geometry'])

linestring = shape(json['features'][0]['geometry'])
geom = Polygon(linestring)

Which leads me to the error

Multi-part geometries do not themselves provide the array interface

If I only use geom = Polygon(shape(json)) I'm getting  Unknown geometry type: featurecollection, which I don't understand – shouldn't it be "FeatureCollection"?

Do I first have to convert to Multipolygon in a Polygon (still in a FeatureCollection) with https://gist.github.com/mapmeld/8742ae89c6d687171d00 ?


texts: How can I find the version of a toolbox used in QGIS

texts: I am using QGIS 2.6.1 Brighton on Xubuntu. The QGIS-install comes with several default toolboxes. Orfeo Toolbox for image analysis is one of them. How can I find out which version of the Orfeo toolbox, or any other toolbox for that matter, is being used with my particular QGIS installation?


texts: PostGIS convert geography(Point,4326) to GPS

texts: I got a database with a column of type geography(Point,4326), the data contained looks like 0101000020E610000092E057F66EF73140930035B56CBD4840. How can I convert it to GPS coordinates?


texts: Getting an ordered list of segments while creating a cross-section through road lines (sidewalks, roadway etc)

texts:
I am trying to get an ordered list of segments of a cross section through road lines. I had a success with Graphical modeler but only in case when road lines were created one by one in one direction (for example, from left to right)

Here is the example of a drawing:

After running a model I am getting a ordered list of segments where I can retrieve start point index, end point index and a length quite easy. The main idea is based on Line Intersections and Points to Path (Processing algorithms).

Here is the package with a project, shp files and a model https://www.dropbox.com/s/4t6fimqfifueww5/transect.zip?dl=0

But when road lines are created without any specific direction (from left to right or right to left), but randomly (for example, first - sidewalks, then centerline, then curbs etc), the Intersection points get ids from these lines and Points to Path produces segments accordingly. So segment 1 is in the middle, then segment 3 is on the right, segment 2 - is on the left etc.
 Here is the package with these shp files https://www.dropbox.com/s/ebwg98isdijasqx/transect2.zip?dl=0

I had a try to use Split with lines the cross-section line, get separate segments (their ids also inherit ids from road layer) and then merge them back. I was hoping to reset the ids, but the merging into one Linestring (all the adjacent vertices were merged too) and a further exploding create segments with the same ids.

Here is a good start for Python, but I am in the very beginning of this process (I was trying to range points of intersection regarding to their distance from the start point of a cross-section line).

Any thoughts and scripts will be appreciated!


texts: Having trouble matching up U.S. Census shapefiles in QGIS

texts:
I'm trying to use the U.S. Census shapefiles to figure out what census blocks make up which state legislative districts. I decided to test the method on a state where I already had this data to see if it was working, and it isn't. QGIS assigned a district to each Census Block, but it was often the wrong one. I'm looking for any help figuring out a good way to get this information.

The process I used was I turned the census blocks into centroids, turned the centroids into stratified random points, joined the original census block shapefile with the random points, then joined this with the state House shapefile. I also tried just doing a spatial join between the block and legislative shapefiles, but had the same problem.


texts: Where is Web Applications in ArcGIS server 10.1?

texts:
I'm new in using ArcGIS server 10.1 . I used ArcGIS server 9.3,9.3.1 and 10 several years.now i want to use ArcGIS server 10.1.When I install ArcGIS server 10.1, I coudln't find Web Applications tab like ArcGIS 10 or later . In ArcGIS server 10 web applications  and web services installed separately but In ArcGIS server 10.1 i just setup "ArcGIS For Server". I could create services but not Web Applications.The questions are that Where is Web Appications in ArcGIS server 10.1 ? How can i create a web application in ArcGIS 10.1 using wizard  ?

Thanks for Advance


texts:
In my industry, KML/KMZ seems to be the standard spatial data format and it is a HASSLE when converting between typical formats.

I have a large number of KML image overlays and am tasked with converting them to TIFF. Typically through QGIS its easy, however I need to automate this given the volume. I've gotten them all written properly, however the image rotations are not preserved.

from gdal import GetDriverByName, OpenEx

tiff_driver = GetDriverByName("GTiff")

in_ds = OpenEx(KML_PATH)

primary_band = in_ds.GetRasterBand(1)

projection = ds.GetProjection()

geo_transform = in_ds.GetGeoTransform()

out_ds = tiff_driver.Create(
     out_path,
     primary_band.XSize,
     primary_band.YSize,
     in_ds.RasterCount,
     primary_band.DataType)

out_ds.SetProjection(projection)
out_ds.SetGeoTransform(geo_transform)

for i in range(1, ds.RasterCount + 1):
     in_band = ds.GetRasterBand(i)
     in_data = in_band.ReadAsArray()
     out_band = out_ds.GetRasterBand(i)
     out_band.WriteArray(in_data)

out_ds.FlushCache()
out_ds = None

The resulting .TIFF files are oriented true north. I can see a "rotation" tag in the KML, however can't find any GDAL options for rotating the new image.

Google Earth screenshot of proper rotation, kml format:

QGIS screenshot of non-rotated geotiff:

I've tried the following:

Replacing the 'rotation' values in the GeoTransform to the KML's rotation in radians
This simple change of the GeoTransform (gt[0], gt1 * cos(rot), gt2 * -sin(rot), gt3, gt[4] * sin(rot), gt[5] * cos(rot))

Neither work properly.

I assume this may involve using the Affine library, however a.rotation(ANGLE) isn't giving any useful results.

Another Update:

Using the Affine library, I can get the proper rotation, however each image seems translated. I'm not manually creating the Affine object, importing from GDAL:

gt = in_ds.GetGeoTransform()
aff_gt = Affine.from_gdal(*gt)

aff_gt *= Affine.rotation(-rotation)

out_ds.SetGeoTransform(aff_gt.to_gdal())

This in effect rotates the original GeoTransform via the Affine library, rotation angle seems legit, however still not giving me what I expected as each image is translated:

Seriously, any direction helps here!


texts: Select By attribute - top 100 values for each district

texts: I am trying to determine the highest 100 values for each district in a shapefile. I have 27 unique districts and 62834 records, but not sure how to proceed. I could "split by attribute" based on the district and select the highest 100 but I would have to repeat this like 30 times for each district.... Is there any easier way or syntax I could use? The value column is number of males in each Small Area Layer in Each district. Using ArcMap 10.3...


texts: Vendor who sells Zip+4 from Semaphore corporation

texts: I need a vendor who sells zip4 CD from Semaphore corporation. I would like to input the address data and zip and then mapped/output that correct address and Zip+4.


texts: How to find number of fields in shapefile using Python?

texts:
I've been looking at the ArcGIS resources for python but I can't seem to figure this one out.
I'm looking for a method to gain access to the number of fields that are within a shapefile / feature class.  Is this information available via Python?


texts: Updating Symbology on Layer using ArcPy?

texts:
I've been trying to update the symbology of a layer on my Python script. The layer that is to be updated shows the percentage of an age range in each region across a county. I have created a test layer (as suggested on GraduatedColorsSymbology ArcPy Website) which copies the symbology of the test layer to the layer I am editing. I have set the symbology of this test layer to use the same dataset as the new layer, but uses a different age range (0 - 16 years).

This worked as expected, so now I want to update the Symbology to use a different age range (65 and Over) using the GraduatedColorsSymbology ArcPy page.

Below is my code:

def Symbology():
    mxd = arcpy.mapping.MapDocument(r"C:\MyMaps\Ages.mxd")
    df = arcpy.mapping.ListDataFrames(mxd)[0]

    updateLayer = arcpy.mapping.ListLayers(mxd, "Age_Range")[0]
    sourceLayer = arcpy.mapping.Layer(r"C:\MyLayers\TEST_LAYER.lyr")

    arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)

    if updateLayer.symbologyType == "GRADUATED_COLORS":
        updateLayer.symbology.valueField = "PCT_65_OVR"
        updateLayer.symbology.numClasses = 5

    mxd.save()

Symbology()

When I execute this code, I get the error ValueError: PCT_65_OVR. I made sure that the attribute name in the code is the same as what is on the symbology layer, so not sure why I am getting this error.


texts: Can't find the proper settings to import Europe map using osm2pgsql?

texts:
I'm trying to deploy a map server, i used http://switch2osm.org/serving-tiles/building-a-tile-server-from-packages/.

My server setup is: Xeon E3-1230v2 Quad-qore, 16GB RAM DDR3, 2x2TB HDD in Software RAID-1 array, running Ubunutu Server 14.04

I use putty to SSH in to it, i don't know how to see why it breaks the import process but i assume the reason is the memory settings. I have dynamic ip at home, so after i set it to import i might lose connection to the server, and when i relog i can't see the import process, but i know it is importing as i check the server load, but it brakes.

Is there a command to show importing process?

I have 7GB of swap memory.

The latest settings i used are:
osm2pgsql --slim  -C 12000 --number-processes 2 europe-latest.osm.pbf

shared_buffers 1024MB

maintenance_work_mem = 1024MB

checkpoint_segments = 100

autovacuum = off (import only setting!)

What should i do?

I would ultimately need the import Europe and North America, also i want to deploy Nominatim 2, for the both regions.

Can my server handle it?

If not what setup am i looking at?

Please be more detailed in your suggestions as i'm not prolific user of Linux.


texts:
I am trying to crop the TIF images based on the shape file. After searching the web, I found the following code and manipulated it.
from osgeo import gdal
import os

input_path = '/.../pics'
output_path = '/.../clipped_pics/'

bandList = [band for band in os.listdir(input_path) if band[-4:]=='.tif']
shp_clip = '/.../Index.shp'

for band in bandList:
    print(output_path + band[:-4]+'_c2'+band[-4:])
    options = gdal.WarpOptions(cutlineDSName=shp_clip,cropToCutline=True)
    outBand = gdal.Warp(srcDSOrSrcDSTab=input_path + band,
                        destNameOrDestDS=output_path + band[:-4]+'_c2'+band[-4:],
                        options=options)
    outBand= None

However, this is the error I get:
TypeError                                 Traceback (most recent call last)
<ipython-input-51-e2d9d6412405> in <module>
      4     outBand = gdal.Warp(srcDSOrSrcDSTab=input_path + band,
      5                         destNameOrDestDS=output_path + band[:-4]+'_c2'+band[-4:],
----> 6                         options=options)
      7     outBand= None

~/opt/anaconda3/lib/python3.7/site-packages/osgeo/gdal.py in Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)
    673
    674     if _is_str_or_unicode(destNameOrDestDS):
--> 675         return wrapper_GDALWarpDestName(destNameOrDestDS, srcDSTab, opts, callback, callback_data)
    676     else:
    677         return wrapper_GDALWarpDestDS(destNameOrDestDS, srcDSTab, opts, callback, callback_data)

~/opt/anaconda3/lib/python3.7/site-packages/osgeo/gdal.py in wrapper_GDALWarpDestName(*args)
   4315 def wrapper_GDALWarpDestName(*args):
   4316     """wrapper_GDALWarpDestName(char const * dest, int object_list_count, GDALWarpAppOptions warpAppOptions, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 4317     return _gdal.wrapper_GDALWarpDestName(*args)
   4318 class GDALVectorTranslateOptions(_object):
   4319     """Proxy of C++ GDALVectorTranslateOptions class."""

TypeError: object of wrong GDALDatasetShadow


texts: Using Plat Plotter - Calculate Curve Table feature given an arc and radius

texts:
Plat Plotter has a Calculate Curve Table feature that generates a sequence of points for use in drawing plot lines. To put another way, it converts a curve into several short, straight segments approximating the curve.

  thence South 13°0'00" West 201.05 feet to the beginning of a tangent
  180.00 foot radius curve, concave Northwesterly; thence Southerly and Westerly along the arc of said curve through an angle of 98° a
  distance of 307.87 feet; thence tangent to said curve North 69°00'00"
  West 255.0 feet to the West line of the Northeast Quarter of the
  Southeast Quarter of said Section 5;

For this example, the straight lines before and after the arc are entered as:

S 13 W 201.05 and N 69 W 255.00

The Calculate Curve Table feature will create a sequence of values that can be placed between the above two entries given a Chord & Radius.

How does one get from the:

  to the beginning of a tangent 180.00 foot radius curve, concave
  Northwesterly; thence Southerly and Westerly along the arc of said
  curve through an angle of 98° a distance of 307.87 feet; thence
  tangent to said curve

to the Chord and Radius values needed by Plat Plotter?

The Chord for Plat Plotter is an angle which will be "S" ddmmss "W" in this case and length fff.ff feet. I assume that the Radius can be taken directly from this legal description, 180.

Chord:  S ddmmss W fff.ff
Radius: 180

So, more specifically, how to calculate ddmmss and fff.ff?


texts: Hide old tiles while loading new ones when zooming in Openlayers

texts:
Standard behavior, while zooming in and out in Openlayers, is that the current tile gets resampled and displayed stretched or compressed until the new tile is available. That's how it looks like while zooming in:

Is it possible to change this behavior on layer level so that the old tile is hidden till the new one is loaded? I would like to apply that to ol.layer.Tile and ol.layer.VectorTile layers.

In my usage case I zoom in quite often when selecting features so tiles from low zoom levels get streched quite heavily. I would prefer to make them invisible and show the background color while loading new tiles.

Note: In Openlayers 2 there was a transitionEffect parameter in OpenLayers.Layer.WMS as shown in this example. 'default' value sounds like what I need, also I don't see any difference to 'resize' value...

As a workaround I could think of listening to tileloadstart and tileloadend events and meanwhile set opacity of the whole layer to 0. But I would prefer to see tiles appearing as the get loaded.


texts: PostGIS: snap latitude/longitude input to nearest OSM road point

texts:
I need a function to return the nearest OSM road point for any latitude/longitude. i.e. users input a latitude/longitude and postgis finds and uses the point on a road that is the nearest to it. Particularly, this needs to be fast.

Here's an extract of my current - slow - solution:

select "geom"
from (select (ST_DumpPoints(geom_vertex)).geom
  from "roads_2po_v") as "vertex"
    order by ST_Distance(geom,
      ST_SetSRID(ST_MakePoint('103.84743100000003','1.2987212'), 4326))
 limit '1')

Unfortunately this order by ST_Distance... limit 1 slows down the entire query about 500x. What would be an improved way to snap to nearest road point in PostGIS using OSM data?

roads_2po_y is a table constructed using osm2po like so:

java -jar osm2po-core-4.8.8-signed.jar prefix=roads cmd=c roads/data.pbf

Structure:

CREATE TABLE roads_2po_v (
  id integer NOT NULL,
  clazz integer,
  osm_id bigint,
  osm_name character varying,
  ref_count integer,
  restrictions character varying
);
SELECT AddGeometryColumn('public', 'roads_2po_v', 'geom_vertex', 4326, 'POINT', 2);

So I'd like to know:

Best/most efficient way to get OSM road points into PostGIS (note: it's for a particular region, not the world)
Efficient query for finding nearest road point when given a latitude and longitude

Any other questions/answers I've found on the site are wrapped in custom requirements. I'm seeking a generic function that returns nearest OSM road point for any given latitude/longitude.


texts: Do Sentinel 2 Level-1C images need atmospherical correction before EVI calculations?

texts:
I would like to use Sentinel-2 Level 1C images to produce EVI maps using the following equation :

EVI=G×(ρnir-ρred)/(ρnir+(C1×ρred-C2×ρblue )+L)

Can I use the Level 1C image as it is, or should I perform atmospheric correction first?

Because EVI is supposed to be an atmospheric resistant index using the coefficent C1 and C2, I'm not too sure if atmospheric correction is necessary or not. My images, from a semi-arid area, are mostly cloud free.


texts:
I am getting the below exception while editing and saving data in spatial table.

Response

<ows:ExceptionReport version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/ows http://x.x.x.x:80/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>{http://x.x.x.x/geoserver/capacity/org}EDIT_TEST is read-only</ows:ExceptionText>
 </ows:Exception>
</ows:ExceptionReport>

I checked the following as suggested in a previous thread:

Primary Key - PK is available in the spatial table and type Number.
User access - On geoserver in stores i have a 'XXXX' user and in oracle i have granted all privileges on the table. But still getting the same issue.

Do we need to do any settings on geoserver?

#

I have noticed that every time new FEATID is generating. Spatial table has primary key and sequence(tableName_colname_sequence). Can any one help me out why the correct featid is not going in the request

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/wfs            http://schemas.opengis.net/wfs/1.0.0/WFS-transaction.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
           <wfs:Delete typeName="capacity:EDIT_TEST" xmlns:capacity="http://x.x.x.x/geoserver/capacity/org">
              <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
                 <ogc:FeatureId fid="EDIT_TEST.fid-270312d3_14ce5aac044_-7ce8"/>
            </ogc:Filter></wfs:Delete>
</wfs:Transaction>


texts: Cannot Save Edits to File Geodatabase Replica in QGIS 2.2

texts:
When trying to edit a File Geodatabase Replica in QGIS 2.2, everything works as expected until the user saves their edits. When pressing save, all edits revert back to their original state with no errors.  If I copy the Feature Class into a new FGDB (thus breaking it out of the replica), I can edit it without any problems.

Anyone have any thoughts on a solution for this?

FGDB Replica created in ArcGIS Desktop 10.2 from a SQL SDE Geodatabase.


texts: Changing coordinate of point using its WKT and PyQGIS

texts:
I need to change the position of a box only if it has a common attribute with a point. In my example if the feature "ID_PARENT" is the same as "CODE".
box = QgsProject.instance().mapLayersByName('box')[0]
point = QgsProject.instance().mapLayersByName('point')[0]

for bfeat in box.getFeatures():
    for pfeat in point.getFeatures():
        if bfeat['ID_PARENT'] == pfeat['CODE']:
            pfeat.geometry() = bfeat.geometry()

Here this code is totally false because I'm assigning a value to a function call but I think it illustrates well what I tried to do which is changing the WKT coordinate.


texts:
Does anybody know of an existing script or a workaround that will convert WFS returned data into something that the ArcGIS JS API can digest and display? Essentially, it would be some sort of script that either returns data that looks like it came from the REST service, or tricks the API into thinking it's an ArcGIS feature layer? I'd prefer not to reinvent the wheel if I don't have to.

Thanks


texts: How to view My Content post on ArcGIS Online?

texts:
I have posted my layers on ArcGIS Online using Esri toolbox and python script.

Now I want to view the layers that are posted.

What is the correct query that will return me Json result of "My Content".

P.S: I can view my layers by logging into arcgis.com and navigating to My content page.


texts: Joining polygons to points when the points fall in the intersection of polygons. QGIS

texts:
I am using QGIS. I am hoping to join a set of geocoded points to polygons in such a way that the new snapefile will be points that have a new field which lists the attributes of the polygon it lies in for each point.  I know this is a simple problem that can be taken care of by the join attributes by location function.

Where the problem comes in for me is that I have many points that lie in the intersection of up multiple polygons.  The outcome I'm looking for is a way to join the points and polygons in such a way that we have points that sequentially list the attributes of each of the multiple polygons it may lie in (with my data the points fall in at most 4 or 5 polygons at once).

Basically the polygons have unique names, I would like to perform some sort of join where I would receive for each point a list of the unique markers of the polygons they fall into.


texts: Calculating proportional area of polygon within another layer's polygon using QGIS?

texts:
I have two polygon layers, model and ntem, whose features overlap and intersect one another. I would like to calculate the proportion of the area of the model polygons which fit within each ntem polygon. Some model features fully fit within ntem features, others overlap multiple ntem features. Both layers are in CRS OSGB:1936.

For example, model feature 50209 falls partly in each of the three ntem features 2202, 2205, 2206 and I'd eventually like to know the proportion of 50209's area in each ntem area.

From initial research, I understand that the Vector > Geoprocessing Tools > Intersect function should do what I want, however I get a completely blank layer of results when I attempt this. (Shapefile is created, empty attribute table, no features!)

Can anyone either point out if I am doing something wrong and suggest an alternative method of achieving the area proportion calculation?


texts: Mask out pixels smaller than 1ha in Google Earth Engine

texts:
I am new here so please forgive any mistakes in this post:
I have a simple script to classify cultivated area vs none cultivated areas using Sentinel-2 in GEE. I want to mask out all objects that are less than 1 ha in size. Below is the script I have put together, but the problem I have is that the .connectedComponents and .connectedPixelCount have a maxSize of 255, which drops very large chunks of the cultivated area's when I finally apply the mask in the end. Can anyone help me tidy this up?

//Add Sentinel imagery from Dry Season 2019, layer stack, composite and clip
var Sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2019-03-20', '2019-08-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
                  .filterBounds(AOI)
                  .select(['B2', 'B3','B4','B8','B11','B12'])
                  .median()
                  .clip(AOI);
//Create an NDVI band and add it to the Sentinel stack
var ndvi = Sentinel2.normalizedDifference(["B8","B4"]);
var with_ndvi = Sentinel2.addBands(ndvi);

// Merge the training points to a single FeatureCollection.
var trainingPoints = Cult.merge(NoCult);

//Obtain signatures for each class from the image
var training = with_ndvi.sampleRegions({collection: trainingPoints, properties: ['Class'],  scale: 10});

//Train the classifier
var classifier = ee.Classifier.cart().train(training, 'Class');

//Run the classification
var classified = with_ndvi.classify(classifier);

//Extract out just the agriculture
var agric = classified.updateMask(classified.eq(1));

// Uniquely label the agric areas.
var objectId = agric.connectedComponents({
  connectedness: ee.Kernel.plus(1),
  maxSize: 255
});

// Compute the number of pixels in each object defined by the "labels" band.
var objectSize = objectId.select('labels')
  .connectedPixelCount({
    maxSize: 255, eightConnected: false
  });

// Get a pixel area image.
var pixelArea = ee.Image.pixelArea();

// Multiply pixel area by the number of pixels in an object to calculate the object area.
var objectArea = objectSize.multiply(pixelArea);

// Mask out objects below 1 hectare
var areaMask = objectArea.gte(10000);

//Create final cultivated area layer with mask applied
var agricFin = agric.updateMask(areaMask);

Map.addLayer(agricFin);
Map.addLayer(agric);

Link to the script is here: https://code.earthengine.google.com/1172c13936c011b39f2c62abf56fecb3


texts: PostGIS extension not available in Postgresql list

texts:
I have installed Postgresql 9.2 from graphical installer in Ubuntu 13.10 and Install PostGIS 2.0.3 through StackBuilder. All things went smooth. But it is not possible to create PostGIS database as there is no extension present in Extension. As far I know two more extension is needed to be there named 'postgis' and 'postgis_topology'. I have tried

CREATE EXTENSION postgis

but it returns there is no postgis.control file within the extension folder. What should I do now?


texts: Creating FeatureCollection from imported assets in Google Earth Engine?

texts:
I have imported several shapefiles as a Fusion Table in Google EE. My goal is to calculate the amount of forest loss (from Hansen et al) and the MODIS VCF (percentage tree cover) for each shapefile, and later do a comparison to see if the two data sets are in agreement. To do so, I want be able to iterate through each shapefile and calculate the amount of pixels of forest loss/VCF per year.

I have one script for the Hansen data set and one for the MODIS VCF, but both are not working ideally for me.

MODIS VCF:

//Import fusion table containing the 8 shapefiles
// Without the ee.filter.eq, I only get a single graph where I don't know for sure from which area it is showing

var fc = ee.FeatureCollection('ft:1Sh6L7D67LeaMWvs9eYnfOQmp-ReR6oQEmWCi2trX')
    .filter(ee.Filter.eq("system:index", '2'))
    ;

var aggregate_array = fc.aggregate_array('name');
print(aggregate_array, fc);

var region = fc;
Map.addLayer(region);

// Zoom to location.
Map.centerObject(region, 7);

// Choose date of interest
var StartDate = '2000-03-05';
var EndDate = '2015-03-07';

// Add the satellite time series
var modisVCF = ee.ImageCollection('MODIS/051/MOD44B')
    .select('Percent_Tree_Cover')
    .filterDate(StartDate,EndDate)
    .filterBounds(region);

print (modisVCF);

// Let's make a chart!
var TimeSeries = ui.Chart.image.seriesByRegion({
    imageCollection: modisVCF,
    regions: region,
    reducer: ee.Reducer.mean(),
    band: 'Percent_Tree_Cover',
    scale: 5000,
    xProperty: 'system:time_start',
    seriesProperty: 'label'
  }).setChartType('ScatterChart')
    .setOptions({
      title: 'Percentage Tree Cover',
      vAxis: {title: 'Tree Cover (%)'},
      lineWidth: 1,
      pointSize: 1,
    });

print(TimeSeries);

My issue with this piece of code is that it only gives me the data I want one location at a time. How can I plot all locations in one graph?

Hansen Global Forest Change:

Which I originally got from here

// Import fusion table
var region = ee.FeatureCollection('ft:1Sh6L7D67LeaMWvs9eYnfOQmp-ReR6oQEmWCi2trX');

// Load Hansen tree cover
var gfcImage = ee.Image('UMD/hansen/global_forest_change_2015');
var lossYear = gfcImage.select('lossyear');

var scale = 100;

//Sequence of years for which we need to extract loss
var years = ee.List.sequence(1, 14); //2001 to 2014

Map.centerObject(region);

var addVar = function(feature) {

  // function to iterate over the sequence of years
  var addVarYear = function(year, feat) {
    // cast var
    year = ee.Number(year).toInt();
    feat = ee.Feature(feat);

    // actual year to write as property
    var actual_year = ee.Number(2000).add(year);
    var year_str = actual_year.format();

    // filter year:
    // 1st: get mask
    var filter = lossYear.eq(year);
    // 2nd: apply mask
    var filtered = gfcImage.updateMask(filter);

    // reduce variables over the feature
    var reduc = filtered.reduceRegion({
      geometry: feature.geometry(),
      reducer: ee.Reducer.sum(),
      scale: scale
    });

    // get results
    var loss = ee.Number(reduc.get('loss'));

    // set names
    var nameloss = ee.String("loss").cat(year_str);

     // alternative 2: always set property
    // set properties to the feature
    return feat.set(nameloss, loss.format());
  };

  // iterate over the sequence
  var newfeat = ee.Feature(years.iterate(addVarYear, feature));

  // return feature with new properties
  return newfeat;
};

var lossstack = region.map(addVar);
print(lossstack);

Map.addLayer(lossstack, {}, "lossstack")

Export.table.toDrive({
  collection: lossstack,
  description: 'ForestLoss',
  folder: "Google EE results",
  fileFormat: "CSV",
  });

In this piece of code, I am able to get all of the forest loss per year per location. However I could use some help in making a more elegant CSV file. Can I choose not to include the geometry somehow?

Here is what the CSV file looks like:


texts: Resolving postgis.control location error on CentOS 7

texts:
Yum installed PostgreSQL 10 and PostGIS 2.0 on a CentOS 7 host. Got this error while creating the postgis extension on a DB:

could not open extension control file
"/usr/pgsql-10/share/extension/postgis.control": No such file or
directory

The postgis.control file on the Centos 7 is actually on /usr/share/pgsql/extension/postgis.control
How should I resolve this? Will creating a soft link be enough?


texts: Indicator Kriging and cokriging using ArcGIS Geostatistical Wizard?

texts:
I have a map of vegetation in 1913 for a particular landscape (approx 35,000ha), of which I have data for just short of 60% of the landscape. The limited visibility is due to the coverage being derived from oblique photos, and large areas are obscured by ridges etc.I need to develop a prediction for what is in the invisible areas. The attached map shows the landscape as Forest (green) , Non-Forest (dark grey) and the missing areas are paler grey (a DEM coverage).

I have settled on using Indicator Kriging (using ArcGIS 10.3 Geostatistical Wizard)  to fill in the missing areas. I have developed a model that predicts forest correctly in 83% of the landscape, and is 86% accurate for non forest (ie. if the real cover is forest, 83% of the time the kriging prediction says it is forest, and 17% of the time it says it is Non Forest).

In my model I have turned Anisotropy to "true", because there is a trend towards more non forest in the SSE corner, and more Forest in the NNW corner. However, I suspect that elevation and insolation are also factors in controlling vegetation.

How can I include Elevation and Insolation (I have the rasters for these) as predictors to see if they improve the model?


texts: Sorting shapefile records on attribute to update new rank field using ArcPy?

texts:
Basically, I am trying to sort the records within the shapefile based on an attribute(eg: pop) and then add a new field that indicates ranking based on the sort. eg: the record with the highest pop would get a rank value of 1, the next would get a value 2.

I can easily sort and add the field in arcgis but the calculating field so that it is ranked based on the order is something i don't know how it works and hence i am thinking may be a python script would be better.


texts:
I want my JSON data with URL can load on the map, and it is a polygon data.
I don't know what the problem with my code.

This is my code:

<link rel="stylesheet" href="https://unpkg.com/leaflet@1.5.1/dist/leaflet.css"
      integrity="sha512-xwE/Az9zrjBIphAcBb3F6JVqxf46+CDLwfLMHloNu6KEQCAWi6HcDUbeOfBIptF7tcCzusKFjFw2yuvEpDL9wQ=="
      crossorigin="" />


texts: How to convert total station data (coordinates) to 3D model?

texts:
i have my coordinates in .csv format and i want to convert it to a usable 3D model.

How do i go about this or what software can i use to achieve this?


texts: How to tell if a country has states and what the state codes are

texts:
With the help of this community I understood how to find country codes and query geoNames data based on that.

What I am stuck on now is understanding how to break things down by state where needed.

For example, take Canada.  There are a few big Provinces like Ontario, Alberta.

I just happen to sort of know that.  But how can I get to know such information about other countries like Mexico or Russia, etc.

Thanks,
Alex


texts: How to rename a geopackage table in Python?

texts:
I'm trying to rename a geopackage table from a QGIS 3.8.0 Python processing script.

But it leaves a broken dataset in 2 parts, where the tabular data table has been renamed but leaves the geometry table not renamed.

Here's the code:

import sqlite3

geopPath = r"C:\GIS\Processing\Tool Processing\NaturalEarth2.gpkg

tab = 'ne_10m_admin_0_countries'

newTab = 'nev_10m_admin_0_countries'

con = sqlite3.connect(geopPath)

cursor = con.cursor()

query01 = "ALTER TABLE '{}' RENAME TO '{}';".format(tab,newTab)

query02 = "UPDATE {} SET table_name = '{}' WHERE table_name = '{}';".format('gpkg_ogr_contents',newTab,tab)

query03 = "UPDATE {} SET table_name = '{}' WHERE table_name = '{}';".format('gpkg_geometry_columns',newTab,tab)

cursor.execute(query01)

cursor.execute(query02)

cursor.execute(query03)

con.commit()

cursor.close()

con.close()

That seems to commit without error. And I suspect query 01 and 02 are doing the same thing.

I don't think editing the gpkg_geometry_columns is needed as it's only a tablename change, the CRS remains the same.

This link is quite useful but can't find the part that helps me.

Which columns in which tables in the underlying geopackage structure do I need to edit?


texts: Retaining symbology used in QGIS when building map in QGIS2web plugin?

texts:
I am a GIS tech for a telco company, I build a lot of maps in QGIS and am always trying to find ways to share the map and its data with my workers.

When i build my map in QGIS and try to open QGIS2web, most of the symbology and just simply the text formating gets all ascew and i dont know why. I've tried using different symbols and changing the formatting but to no avail.

Here is the map in QGIS:

And here's the output in QGIS2Web:

the 1 and 2 you see are telephone poles I want people to be able to click and get the relevant information, but as it stands the icons or text i use simply doesn't show up.

Any advice you have with working with this tool , or offering an alternative is perfectly acceptable. the only clincher is that it needs to be able to work locally.


texts:
Is there a tool that will pretty-print well-known-text strings? One that would accept for instance

MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 45 20, 30 5, 20 35)))

and output something like

MULTIPOLYGON (
  (
    (
      40 40,
      20 45,
      45 30,
      40 40
    ),
  ),
  (
    (
      20 35,
      45 20,
      30  5,
      20 35
    )
  )

Ideally this would be an online tool that runs in a browser (along the lines of the many JSON formatters that exist), but an offline tool would still be useful.

There are a number of online tools for viewing and editing WKT-format data, but none that seem to offer any support for reformatting. The OpenLayers vector-formats example app has a checkbox tantalizingly labelled "Pretty print", but selecting it seems not to do anything when WKT is set as the output format.


texts:
I am writing a script to use an updatecursor.  The feature class is locked so I am versioning the database,  creating a feature layer, and reconciling versions.  I have taken the skeleton for my script from this thread. Debugging RuntimeError: workspace already in transaction mode from arcpy.da.UpdateCursor and ArcSDE feature classes?

The main difference is that I already had an SDE connection file in my directory. Should this be re-created each time the script is executed to avoid this error?

__author__ = 'Administrator'

from arcpy import *
import GetSRConfig
import os
import time

def GetDateTimeString(n = None):
        """ format a datetime to string """
        if(n==None):
            s = time.strftime("%Y%m%d%H%M%S", time.localtime())
        else:
            s = time.strftime("%Y%m%d%H%M%S", time.localtime())
            if((isNumeric(n)==True) and ((n>4) and (n<14))):
               s = s[0:n]
            else:
               s = s[0:14]
        return s

def isNumeric(s):
    b = True
    try:
        i = float(s)
    except:    # not numericelse:    # numeric
        b= False
    return b

#Locals
P6featureName = GetSRConfig.SDELayer
Parent = "dbo.DEFAULT"
version = "myVersion" + GetDateTimeString(12)
featureLyr = "lyr" + GetDateTimeString(12)

print version
# Server = ***
# Service = ***
user ="dbo"
# Pass = ***
# SDE = 'E:\C_Drive_files\Administrator'\311Request\data\ServiceRequest.sde'
temploc = ""
fil = "SDETempConn"

env.overwriteOutput = True

#Create Version
print "Creating version"
arcpy.CreateVersion_management (GetSRConfig.SDEConnFile, Parent, version, "PUBLIC")
VersionName = user.upper() + "." + version

#Create new connection
workspace = os.path.join (temploc, fil + ".sde")
print workspace
#Layers
P6feature = os.path.join (workspace, P6featureName)
arcpy.MakeFeatureLayer_management (GetSRConfig.SDELayer, featureLyr)

#Start editing
print "Initiating editing"
edit = arcpy.da.Editor (GetSRConfig.SDEConnFile)
edit.startEditing ()
edit.startOperation()

#Test Cursor
print "Testing cursor"
P6Cursor = arcpy.da.UpdateCursor (featureLyr, ["VehicleNam"])
for row in P6Cursor:
    print row[0]
del row
del P6Cursor

#Stop/save edits
edit.stopOperation()
print "Stopping editing"
edit.stopEditing("True")

#Switch to version
print "Switching version"
arcpy.ChangeVersion_management(featureLyr, "TRANSACTIONAL", Parent)

#Reconcile and post
print "Reconciling and posting"
arcpy.ReconcileVersions_management (GetSRConfig.SDEConnFile, "", Parent, VersionName, with_post = "POST", with_delete = "DELETE_VERSION")


texts: How to visually join river to lake in QGIS?

texts:
I have two shapefiles, one is for landuse (which contains forests, and also lakes), and another one for rivers. When I set the both water types to have outlines, they are not joined nicely - they overlap, because river end has an outline too. Is there a method to join them "nicely", like you do with roads (and enabling levels)?

Example here:


texts: Why are Bounding Coordinates not as expected?

texts:
I have a shapefile in GCS_WGS_1984 projection, but the bounding coordinates/Feature extent shows the wrong coordinates in the screenshot below.

.

I tried to change the projection, by deleting the *.prj file and defining a new projection but it still did not get a good result.

Can anyone explain what I may be doing wrong?


texts: ArcGIS javascript API setDefinitionExpression returns wrong result set

texts:
New to GIS and ArcGIS.

I am trying to use the javascript API setDefinitionExpression to limit the output of a FeatureLayer on the map; however, the result set returned by the function is incorrect. The function returns only one result set but when I try to run the same query on the either the ArcGIS server or the database server, there are 6 results.

Anyone know why this is the case?

I am using SQL Server 2012 and following is the code I am using:

var layer1 = new FeatureLayer(path + "/arcgis/rest/services/MammalDistributionMapOIMap/MapServer/0", {
    outFields: ["*"],
    opacity: 0.5
});
var speciesID = document.getElementById("speciesID").value;
stringquery = "species_id = " + speciesID;
layer1.setDefinitionExpression(stringquery);


texts: How to create 2D or 3D cross sections from geological borehole data in QGIS

texts:
Is it possible to load borehole csv data that contains x, y, z coordinates into QGIS and create cross sections of the strata? It is data along a rail line so a DEM is not suitable. Just a 2D cross section would be fine.

I have looked through the similar question threads but the problem was never properly answered.


texts: Making Python script tool run immediately upon opening?

texts: I have an ArcGIS 10.1 script tool that creates an ArcSDE connection folder and copies data from the server to a new folder on the users computer. The tool has no parameters runs fine but when you open the tool it gives the default "this tool has no parameters" interface the user then presses ok and the tool runs. Is there anyway(perhaps with toolvalidator) to just run the script when the tool is opened without this pop-up?


texts: How to avoid empty popups when GetFeatureInfo request is empty

texts:
I set up a WMS with QGIS SERVER. I am querying it via a web app based on OpenLayers.

Everytime I click on my map, I send a GetFeatureInfo to a specific layer passing the coordinates where I clicked, and the response is rendered in a popup. So far so good.

The problem is if I click where there's no layer's feature in the map, I'd like no popup to show up, but instead I get an "empty-table-response" from the WMS which is in turn rendered in my popup, and this is annoying.

Here is an example of a request with some interesting data:

https://www.wondermap.it/cgi-bin/qgis_mapserv.fcgi?map=/home/ubuntu/qgis/projects/Demo_sci_WMS/demo_sci.qgs&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=piste_sci&LAYERS=piste_sci&INFO_FORMAT=text%2Fhtml&I=50&J=50&CRS=EPSG%3A3857&STYLES=&WIDTH=101&HEIGHT=101&BBOX=1164063.3256063634%2C5842307.561011146%2C1165993.360570564%2C5844237.595975346

And here is one with empty data which I don't want to display in the popup:

https://www.wondermap.it/cgi-bin/qgis_mapserv.fcgi?map=/home/ubuntu/qgis/projects/Demo_sci_WMS/demo_sci.qgs&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=true&QUERY_LAYERS=piste_sci&LAYERS=piste_sci&INFO_FORMAT=text/html&I=50&J=50&CRS=EPSG:3857&STYLES=&WIDTH=101&HEIGHT=101&BBOX=1161292.4833310258,5842938.166494499,1163222.5182952264,5844868.201458699

This is how I am requesting the data in OL (only the relevant part):

function getInfoUrl(evt, view, lyr) {
    var coordinate = evt.coordinate;
    var url = CORS_PREFIX;
    var viewResolution = /** @type {number} */ (view.getResolution());
    if (lyr.get('type') != 'base') {
        var visible = lyr.getVisible();
        if (visible) {
            url += lyr.get('source').getGetFeatureInfoUrl(evt.coordinate, viewResolution, 'EPSG:3857', {
                'INFO_FORMAT': 'text/xml'
            });
            return url;
        }
    }
}

map.on('singleclick', function(evt) {
        // display info in fancyAlert
        for (let i = 0; i < ol_layers.length; ++i) {
            let layer = ol_layers[i];
            let text = getInfoUrl(evt, view, layer);
            if (text) {

                var xhttp = new XMLHttpRequest();
                xhttp.onreadystatechange = function() {
                    if (this.readyState == 4 && this.status == 200) {
                        fancyAlert(this.responseText, 'info', 'Layer Info');
                    }
                };
                xhttp.open("GET", text, true);
                xhttp.send();
            }
        }
    })

How should I avoid the second response to show up? I thought a request to an empty feature area should trigger no response or at least a WMS exception. Maybe this is not a standard but it's what I was used to.


texts: Upgraded ESRI .NET ADF application (9.3 to 10.1) intermittently missing service issue

texts:
I also asked this question on the ESRI forum but I'm not sure if anyone will see it through the tumbleweeds over there.

I had an application based on the .NET Web ADF originally built in 9.3. The server is being upgraded to 10.1 so the application needs to move with it (re-architecting is not an option).

The automated upgrade utility for moving ADF applications to 10.1 was not the smoothest process and I have since fixed a couple of very obvious bugs in the ESRI JavaScript. However there is one problem I can't figure out:

Intermittently, after a full page reload, the map fails to render. The application calls ESRI.ArcGIS.ADF.Web.UI.WebControls.MapHandler.ashx requesting a URL for the new map image. Sometimes that handler returns

{"code":400,"type":"error","errorMessage":"Resource <service name> not found"}

The application then requests an image from ESRI.ArcGIS.ADF.Web.MimeImage.ashx passing an ImgID parameter that seems to be randomly made-up - I can't see it in any previous requests and it obviously isn't provided by MapHandler.ashx. The server then responds with a 404. If I reload the application again the same requests are sent with exactly the same parameters (in the example I saw the post-reload request even passed the same 't' parameter which I assume should represent time) and a successful response is returned.

Can anyone tell me why my service is intermittently reported as unavailable? ArcGIS Server is not logging any errors or warnings. Based on this and how quickly the request is re-issued and successfully answered I don't think the service is actually having trouble. Comparing a working and a failing request in Fiddler gives me no clues as the same headers and cookies seem to be sent.

EDIT
The only reason I can imagine is if the server control 'MapResourceManager' or something associated with it isn't fully initialised by the time the first request is sent, but as such timing issues are the responsibility of the framework I would be surprised if this was the problem. However I've seen widespread use of setTimeout in the ADF's JavaScript and I'm really hoping this issue isn't the result of some flawed assumption about timing.


texts: Identifying Europe's major road exits using QGIS

texts:
Dataset from EuroGeographic and Global Roads Open Access Data Set (gROADS)  have been downloaded.
The goal is to identify the exits (geocoordinate points) of highways and primary roads (gROADS refers to them as FClass     Functional Class            1=Highway, 2=Primary).  The shape data is seen and visible with software tool QGIS, but the roads are but collections of lines.
The problem is the data set is not obvious to parse.  The EGM 2021 Dataspecification has a Data Dictionary; under the transportation section data tables a 'Feature class name' EXITC  with code  AQ090 which are not found or, if present, sometimes of limited data. So for example, country code IT has no Exitc dbf file.  Yet Google and Openstreetmap have these elements.
Consulting the database files (EuroGeographic) for ROADL which are lines, entries have names that are repeated, thus I assume that these are segments, and thus the end-points of the segments would correspond to entry/exit points.
Yet I cannot fathom how these can be extracted from this data example.
OBJECTID    FCsubtype   inspireId   beginLifes  F_CODE  ICC COR EXS LLE LTN MED NAMN1   NAMN2   NAMA1   NAMA2   NLN1    NLN2    RST RSU RTE RTN RTT TEN TOL TUC SHAPE_Leng  geometry
3   1   ROADL1664975    20180614    AP030   IT  1   28  2   4   1   Autostrada dei Laghi    N_A Autostrada dei Laghi    N_A ITA N_A 1   1   E35 A9  16  1   2   7   0.00262150576963    LineString

How can these exits be obtained?


texts: OpenLayers Clustering Adds Features to Selected Features Array on Zoom

texts:
This is my first posted question despite using StackExchange for a few years now. I've been stuck for quite some time trying to figure out this issue with OpenLayers. The problem is as follows:

I have a bunch of features plotted on a map and am able to select features in various ways (e.g. drawing a box around them, clicking on individual features, entering a zip code that selects all features within that zip code, etc.). This all works fine. The issue is that I have a clustering strategy on the features layer and when I zoom out, currently selected features get clustered together with unselected features and those unselected features will become selected when I zoom back in. What I'm trying to accomplish is to NOT have the previously unselected features become selected after clustering.

Here is a link to a visual hosted on my website (http://korynorthrop.com/uploads/Visual.jpg) I could post some code as well, but thought I'd simplify this question first.


texts: How to properly combine and export (as raster) vector lines on top of raster

texts: I have a georeferenced TIFF with building plans on it. I then laid topo lines on top of it. it looks right in QGis. I now want to export this as an image (i.e. with the topo lines overlaid). How do I do this? Thanks


texts: Clip/Intersect in QGIS 2.8.2 resulting in empty attribute table

texts:
I'm trying to clip two layers together (a layer showing the area of a city and a layer showing a political boundary running through that city) in order to show the location of the political boundary within the city limits only. I've attempted to achieve this using both the 'clip' and 'intersect' tools, and have gotten the same result with both--the layer is listed as existing but is not visible on my map, and the attribute table for the resulting layer is empty. I have tried turning off 'on-the-fly-reprojection,' I've made sure the layers are using the same projection, and I'm definitely using the correct input and clip/intersect layers. The result is always an empty layer.

Based on reading the answers to similar questions, I think that running a spatial query might get me the result I need.


texts:
I have managed to open and intersect a gridded .nc file with spatial polygons for river basins. Unfortunately, I cannot find which variable indicates the month at when the values were measured. The file should contain monthly irrigated area data.

Here is my code to open the file

irrig_path <- "/Volumes/Transcend/Uni/ETH/International Political Economy/data/Huang Water Use/irrigation water use/"
irrig_name <- "cons_irr_h08.nc"
irrig_file <- paste(irrig_path, irrig_name, sep="")
irrig_r <- raster(irrig_file, varname="cons_irr")
irrig_r

Here is a link to the file

https://drive.google.com/file/d/1Ry68E_M3e0twpBn8V8NdwSsHKkERE3AR/view?usp=sharing

and a link to the article that published the data. I could not find any indication where the date is stored.

https://www.hydrol-earth-syst-sci.net/22/2117/2018/hess-22-2117-2018.pdf

This is an image of the content of the raster data once loaded into R


texts: How to calculate point(or XY) for Start point of Slope and Endpoint of Slope along the Section line on DTM. for instance


texts: SCP - Keyerror when merging ROIs into one class

texts:
I want to merge the spectral signatures from the classes in the ROI list into one class. E.g. I have 4 ROIs representing water and I want to merge these 4 water ROIs into only one water class to continue my classification with SCP.
But regardless the ROIs or the correspondig type I get this message everytime:

QGIS Version: 3.16.4-Hannover
Ubuntu 18.04


texts: How to install stable GDAL on Ubuntu 11.10?

texts:
I am trying to install the stable ubuntu package from here https://launchpad.net/~ubuntugis/+archive/ppa/+packages on my ubuntu 11.10 server

I have added the ubuntugis stable repository to server's list of software sources with these commands

sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ppa

But, apt-get install libgdal1-1.8.0 command fails to find the package. Is the ubuntu 11.10 server not supported? I am newbie with this, so i am not sure if the missing oneric (ubuntu 11.10) entries is a sign that it is supported.


texts: How to add href to popup (infoWindow) in ArcGIS API?

texts:
I'm a newbie in programming. I'm trying to make a href link inside of arcGIS js infoWindow (popup). I don't know how to make this link work. The code that I'm trying is shown below.
Could you please suggest a solution?

    map.on("click", addPoint);
    function addPoint(evt) {
        var latitude = evt.mapPoint.getLatitude();
        var longitude = evt.mapPoint.getLongitude();
        var link = "Google Office"
        map.infoWindow.setTitle("Coordinates");
        map.infoWindow.setContent(
          "lat: " + latitude.toFixed(7) + "
 lon: " + longitude.toFixed(7) + "
 Link: " + link
        );
        map.infoWindow.show(evt.mapPoint, map.getInfoWindowAnchor(evt.screenPoint));
      };


texts:
I'm new to GRASS and GIS environments in general, and I recently faced a problem when trying to get the profile plot from a DTM. When I draw the transect, instead of the plot, I get the following message:

Traceback (most recent call last):
  File "C:\Program Files (x86)\GRASS
6.4.2\etc\wxpython\gui_modules\profile.py", line 468, in
OnCreateProfile

self.DrawPlot()
  File "C:\Program Files (x86)\GRASS
6.4.2\etc\wxpython\gui_modules\profile.py", line 490, in
DrawPlot

marker=self.properties['marker']['type'])
  File "C:\Program Files (x86)\GRASS 6.4.2\Python27\lib
\site-packages\wx-2.8-msw-unicode\wx\lib\plot.py", line 342,
in __init__

PolyPoints.__init__(self, points, attr)
  File "C:\Program Files (x86)\GRASS 6.4.2\Python27\lib
\site-packages\wx-2.8-msw-unicode\wx\lib\plot.py", line 145,
in __init__

self._points =
_Numeric.array(points).astype(_Numeric.Float64)
ValueError
:
could not convert string to float: *

Does anyone out there have an idea about my problem?


texts: Using row.setValue with arcpy.da.UpdateCursor gives 'int' object has no attribute 'setValue'

texts:
I am adding a new field to my attribute table, and trying to populate it with a a section that is already in the attribute table but converted to km instead of inches. My code does not work with either the index [1] or [0] giving me the error

  list index out of range

and

  'int' object has no attribute 'setValue'

Here is my code:

target= ["C", "RW", "LW"] #positional query
height = "height"
weight= "weight"

arcpy.MakeFeatureLayer_management("Roster", "Center", "position = '"+target[0]+"'")
arcpy.AddField_management("Center", "heightkm", "FLOAT")
with arcpy.da.UpdateCursor("Center", "height_km") as cursor:
    for row in cursor:
    row.setValue("height_km", row.getValue("height") * 2.54)
    cursor.updateRow(row)

I am unsure what I am doing wrong.


texts: Selecting and extracting data from TIFF files

texts:
I have some tiff files (2 for example) and a shapefile.
The real format of the tiff files is geoTiff.

I want to get the points inside the shape boundary. Each point should have 4 parameters: geom(coordinates), red, green, infrared (we have tiff files with such info).

At the end of the operation I want to get a vectorized shape with points.

Could you please give me some advice on how can I do that programmatically.
Are there any libs/tools?
It is all in the context of geoserver/postgresql/geomajas framework.

Upd: now I have more information. My tiffs will be one of the Landsat combination: http://gif.berkeley.edu/documents/Landsat%20Band%20Information.pdf

The main question now is: to get RGB channels info in exact point, should I parse image points on RGB,or I can get this info from geoTiffMetadata (there is http://download.deegree.org/deegree3/nightly/core/javadoc/org/deegree/coverage/raster/io/imageio/GeoTiffIIOMetadataAdapter.html ). I didn`t find really useful info while debugging the code and calling the methods of GeoTiffIIOMetadataAdapter.


texts: Does ArcGIS Desktop 10.7.1 work with Python 3.9?

texts:
I wish to install ArcGIS 10.7.1 on my laptop running Windows 10.  It appears to install Python 2.7 as part of the process. However, I have Python 3.9 installed for prior course work.
If 10.7.1 requires 2.7 of Python, can both versions of Python coexist or will it mess up all my environment variables?
I also run Jupyter-Lab the latest version.


texts:
I have a locally hosted Nominatim service. I'm trying to reverse geocode a lat, lon to a place hierarchy.

E.g.
https://nominatim.openstreetmap.org/search.php?q=30.284780%2C+-97.738831&viewbox=-126.57%2C51.68%2C126.57%2C-51.68

which returns https://nominatim.openstreetmap.org/details.php?place_id=93023332

How do I travel up the place hierarchy from http://www.openstreetmap.org/way/150732423 to http://www.openstreetmap.org/relation/1701848#map=12/30.3361/-97.7536 ?


texts: How to convert 6 & 7 digit coordinates to degrees latitude & longitude

texts:
I'm a bit new to the whole GIS world, and though I am making great headway drawing polygons and such in Bing Maps, I am currently stuck. In our SQL spatial data, I came across MULTIPOLYGON data, where all the coordinate values are very large numbers (6 digits and 7 digits for what should be latitude and longitude, before the decimal). I found an SQL method to dismantle the multi-polygon into several polygons, but those large values do not work in Bing Maps nor Google Maps. So I've been searching for an algorithm to convert those large numbers down to what I have perceived as regular latitude/longitude numbers, such as 45.19167 -93.58328 (which does work on online maps).

With that said, could somebody please either show me or point me in a great direction to learning how to convert 631456.83, 3816200.79 to something I can use online. I know the coordinates fall somewhere in Arizona.


texts: Hardware and software for capturing bathymetric lake elevations?

texts:
What hardware and software would you use for capturing bathymetric elevations (depth) in inland waterbodies (lakes and rivers)?

Your device and software budget is $1,000. The budget is a soft target, to get setup, and doesn't include actual data capture. A 17' open boat with outboard motor and a ruggedized laptop running Windows XP Pro and Arcmap (Arcview license) is at your disposal, if you need them. The data captured should be suitable for building a raster elevation model and bathymetric contours. Depths are not expected to exceed 200m. Desired accuracy should be comparable to or better than a standard consumer handheld GPS (Garmin GPSMap, etc.).

(Actually the budget is a sham. If you've got a good solution for $5k or $50k I'd like to know about it too.)


texts:
I've been looking at making an azimuthal equidistant map. There are several web sites that offer to make one, based on my location, but they all have the same problem. They assume my antenna rotates 360 degrees - so they map the entire world, staarting from my location, and finishing at 20,000km  (1/2 earth circumference).

However, my antenna setup only covers 180 degrees.
As such, I'm looking for a "half circle" shaped map, that ranges to 40,000 km.

Does such a map exist? Is "azimuthal equidistant" simply the wrong search term, leading me astray?


texts: Install PostGIS from source for specific PostgreSQL version

texts:
I have two versions of PostgreSQL on my machine (Ubuntu 12.04) one is 9.1 and the other one is 9.4.

I want to install from source (I have my reasons) postgis 2.1.2 over PostgreSQL 9.4, but when I ran configure it finds only 9.1.

How to specify 9.4 version during compilation process?


texts: Projecting selection from large data set using ArcPy?

texts:
I have a large data set X(x1-xn) from which I must select consecutively the attributes A(a1-an) from X.fieldsSRC or X.fieldDST and project the selections.

I have a functioning arcpy script that does the following:

<inside loop incrementing n>
#proj_sr is defined here but omitted
where_clause = fieldSRC = n OR fieldDST = n #simplified for webparser

#create a new feature class based on selection
arcpy.Select_analysis('input', 'temp1', where_clause)
#project selection
arcpy.Project_management('temp1', 'temp2', proj_sr)

#PERFORM OPERATION ON PROJECTED SELECTION

#append selection, and project to destination SR
arcpy.Append_management('temp2', 'final_output')

#we must delete rows from source files otherwise we will treat them each twice
#create layer file in order to select and delete specific rows
arcpy.MakeFeatureLayer_management('input', 'layerfile', where_clause)
#delete rows selected above in layer file
arcpy.DeleteFeatures_management('layerfile')
<end of loop>

This works. In my mind there should be a manner to do the following which would be more efficient.

<inside loop incrementing n>
#proj_sr is defined here but omitted
where_clause = fieldSRC = n OR fieldDST = n #simplified for webparser

#create layer file in order to select and delete specific rows
arcpy.MakeFeatureLayer_management('input', 'layerfile', where_clause)
#project ACTUAL selection rather than a complete shapefile as was done earlier
arcpy.Project_management('layerfile', 'temp2', proj_sr)

#PERFORM OPERATION ON PROJECTED SELECTION

#delete rows selected above in layer file
arcpy.DeleteFeatures_management('layerfile')
<end of loop>

Obviously the arcpy.Project_management() function does not operate in the second code example.
It's a rather large data set ('input') and recreating the feature class based on the selection each time does not seem efficient/necessary. It works but it is not fast.

Does anyone have a suggestion?


texts: Layer displaying correctly with OpenLayers and GeoServer as WMS but giving CORS error when trying to display the same layer as WFS

texts:
I have the following map which works perfectly:

import 'ol/ol.css';
import {Map, View} from 'ol';
import TileLayer from 'ol/layer/Tile';
import TileWMS from 'ol/source/TileWMS';
import BingMaps from 'ol/source/BingMaps'
import VectorLayer from 'ol/layer/Vector';
import VectorSource from 'ol/source/Vector';
import GeoJSON from 'ol/format/GeoJSON';

var baseLayer = new TileLayer({
source: new BingMaps({
   key: 'mykey',
   imagerySet:'AerialWithLabels'}),
});

var map = new Map({
  layers: [baseLayer],
  target: 'map',
  view: new View({
    zoom: 11,
  }),
});

When I add a WMS layer to this map it displays correctly:

var layer = new TileLayer({
  source: new TileWMS({
    url: 'http://myip:8080/geoserver/MYWORKSPACE/wms',
    params: {'LAYERS': 'MYWORKSPACE:MYLAYER', 'TILED': true},
    serverType: 'geoserver',
  }),
});

map.addLayer(layer)

But when I add the same layer but as WFS:

var layer2 = new VectorLayer({
  source: new VectorSource({
    format: new GeoJSON(),
    url: function(extent) {
      return  'http://myip:8080/geoserver/MYWORKSPACE/wfs?' +
              'service=WFS&' +
              'version=1.0.0&request=GetFeature&typename=MYLAYER&'+
              'outputFormat=application/json';
    }
  })
});

map.addLayer(layer2)

It shows nothing and gives me the following error:

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource at
  http://myip:8080/geoserver/MYWORKSPACE/wfs?service=W…GetFeature&typename=MYLAYER&outputFormat=application/json.
  (Reason: CORS header ‘Access-Control-Allow-Origin’ missing).

I know I could enable CORS on my server and maybe it would work but I don't understand why displaying the layer as WMS doesn't give me this CORS error.

I wouldn't like to enable CORS on my server because I already tried that last day without success, so If there was a way to show this layer as a WFS without touching the server would be perfect.


texts: Importing processing from QGIS into IDLE

texts:
I have been trying to make a program that creates a virtual raster based on the files that are input. Whenever I run my code it says.

Traceback (most recent call last):   File
"C:/Users/Zach/Desktop/Python_testing/Geo_test_1.py", line 1, in

import processing ModuleNotFoundError: No module named 'processing'

My code is:
import processing

file_input=['C:/Users/Zach/Desktop/Qgis/Missisipii_Delta/LC08_L1TP_021040_20201011_20201016_01_T1_B1.TIF'.... there are 7 more bands]

file_output = 'C:/Users/Zach/Desktop/Qgis/proc_build.TIF'
processing.run("gdal:buildvirtualraster",\
{'INPUT':file_input,\
'RESOLUTION':0,'SEPARATE':True,\
'PROJ_DIFFERENCE':False,\
'ADD_ALPHA':False,\
'ASSIGN_CRS':None,\
'RESAMPLING':0,\
'SRC_NODATA':'',\
'EXTRA':'',\
'OUTPUT':file_output})

I am not sure if this is the problem, but I installed GDAL through Anaconda, and Anaconda overall has seemed to bring a few complications. I do not know where to find processing in my directory, nor how to make it so that it works in IDLE. This code works fine in the QGIS Python Console. So, how to import processing from QGIS into IDLE?


texts: Creating separate viewsheds for individual points using ArcGIS ModelBuilder?

texts:
I have a feature class with 3000 points. My goal is to have a viewshed, at each point, clipped to 50m and 500m buffers. I want to be able to look at the percent visible/not visible for each point. I found the following link explaining how to use an iterator, but it's not super explicit and my attempt to copy it is not giving me the results I want. Here is the link: https://geonet.esri.com/thread/31124

The first answer has an attached model I tried to copy. I have attached a png of my model.

My results are one raster file with all of the buffers (so it kind of worked), it doesn't look like it did it for each individual point because only one large cluster of maybe a few hundred points have both visible and non visible values.

I'm not very familiar with iterators and only moderately familiar with ModelBuilder, is there something obvious you can see that I'm missing?


texts: Specify geoprocessing service connection at publish

texts:
I have a geoprocessing service published from a tool written in Python and contained inside a Python toolbox.

I need to be able to publish this service in multiple environments (production, staging, testing, development) i.e. using different database backends.

I do not want to recreate the service definition package every time I need to move to a different environment.

In particular, I do not have direct access to production and do not want the people deploying the app to need to fire up ArcMap or ArcCatalog and risk some parameter not being set correctly.

Instead, I want to hand them a package and say "Publish this using the connection you created".

How can I accomplish that?

Edit:

Thanks to another question, I've run across this piece of documentation here:

  If you are currently disconnected from the server or do not have
  access to a server connection, you can configure your service
  definition with no available connection to ArcGIS Server. Selecting
  this option creates a service definition file that will have to be
  configured to work with a server connection at the time of publishing.

So it seems to me that ESRI says this is possible, unless I've misinterpreted something.


texts: Changing local coordinate system to WGS84 in QGIS print composer?

texts:
I have composed a map where shapefiles, points etc. all are in local coordinate system (l-est).

Is it possible to convert it to common WGS84 system?

So that in printed map there are coordinates 23.906250;59.432234 (not local 6506050; 4534232)?

I've tried to change it under items properties -> the CRS selector, but no luck.

At the moment in "QGIS print composer" I have this picture:


texts: Export Raster from ArcMap gdb raster catalog

texts:
I received a file geodatabase raster catalog and would like to extract the rasters as GeoTIFF.

I tried to read them in FME and save them as raster but it failed and I also tried to export them in ArcMap (select all polygons in the layer and -->data-->export but this seems to export all rasters as a single raster.

How can I batch export all raster images which appear in ArcMap when I zoom in so I get a folder of GeoTIFFs from this raster catalog?


texts: I have a 12 mb KML file that I can easily import into QGIS.  The trouble is when I save it as a shapefile the data in the description field, which is html, becomes truncated.  Is there a way to avoid this?


texts: Adding M values to point feature attribute table in Arcmap

texts:
I want to create a point feature with XYZ and M values. The points are 10m spaced points located on a line. I need this to create cross sections in Grapher. My problem is I cannot seem to get the M values by "add geometry attributes" tool. I first created points along a line(every 10m). Then I made vertices along the line and converted the vertices to points, but the distance along the line just does not get carried over.

ArcMap 10.4.1


texts:
I've created a geopackage through gdal_translate (overviews with gdaladdo)

In geoserver 2.13/2.15 I can create the geopackage store with no problems, but then geoseverver won't fetch the tables/layers from the geopackage.

This is the geopackage gdalinfo

 Driver: GPKG/GeoPackage
Files: ortos2017.gpkg
Size is 133333, 133333
Coordinate System is:
PROJCS["ETRS89 / Portugal TM06",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",39.66825833333333],
    PARAMETER["central_meridian",-8.133108333333334],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","3763"]]
Origin = (-33600.000000000000000,83000.000000000000000)
Pixel Size = (0.120000000000000,-0.120000000000000)
Metadata:
  AREA_OR_POINT=Area
  IDENTIFIER=ortos2017
  ZOOM_LEVEL=10
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -33600.000,   83000.000) (  8d31'44.35"W, 40d24'54.32"N)
Lower Left  (  -33600.000,   67000.040) (  8d31'41.32"W, 40d16'15.61"N)
Upper Right (  -17600.040,   83000.000) (  8d20'25.71"W, 40d24'56.09"N)
Lower Right (  -17600.040,   67000.040) (  8d20'24.12"W, 40d16'17.37"N)
Center      (  -25600.020,   75000.020) (  8d26' 3.88"W, 40d20'35.99"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 66667x66667, 33333x33333, 16667x16667, 8333x8333, 4167x4167, 2083x2083, 1042x1042, 521x521

What am I missing? Do I have to add other attributes to the gpkg?


texts: Deploy SD file to production environment

texts:
I have created map service using an oracle 11g geodatabase.
I have created the SD file which knows the exact location of my table "myPlaces". (IP of my XE, credential etc)
Now I have to send the SD file to my client so as him to publish it.
Supposing that he has created his geodatabase (in the production environment) and he has created his own user/ credentials.
And Suppose that he has registered  his geodatabase with his server.

How Will the service undestand that the it will retrieve the data from the table "mylocations" of his geodatabase, since I cannot specify it in some way


texts: Android ArcGIS API, how to loop into every Feature in Geodatabase

texts:
I have an offline geodatabase in android ArcGIS which I have downloaded through gdbAsyncTask.

how do I loop into every feature in that geodatabase?

for example I could have collect all the IDs from the geodatabase and then loop into each in a for loop but this is not possible.

new FeatureLayer(localGdb.getGeodatabaseFeatureTableByLayerId(id)).getFeatureIDs(x,y,tolerance);

getFeatureIDs method does not work since I guess I have to also add the layer in my map. But I want to do that without having to add it to my map. Since I need to get the geometries only.

After the point: localGdb.getGeodatabaseFeatureTableByLayerId(id), how do I guery for all the features? or query for all the feature IDs, so then I can reach every feaure?


texts: GeoServer Import of MBTiles file Results in Garbled Map

texts:
I created a mbtiles file using MapTiler 9.1, and imported it to GeoServer (no errors). When I previewed the layer, it came out all garbled up (see screenshot below). I could see this being an issue with the data import settings, but have no idea where to begin to address this. I don't think the images themselves are corrupted, as they preview fine in MapTiler.

Oddly enough some tiles seem to render fine at different zoom levels, some garbled, and others not at all.

Update: after playing around with Layer Publishing settings, it appears to do with the Interpolation Method. Unfortunately any setting I chose resulted in better or worse images, leading me to believe that the export from MapTiler used some awkward interpolation settings. I am re-exporting now using Nearest Neighbor settings, which is also present in GeoServer.

Update 2: the new export didn't change anything, the behavior remained the same in GeoServer. I have the default raster style applied to this layer. Seems this might be an issue/incompatibility with the rendering engine?


texts: Setting symbol diameter based on specific attribute parameter in ArcGIS Desktop?

texts:
I am working with a shapefile of a forest, and I want the symbol of each of the trees to reflect their exact diameter.

I have tried to sort this out through the Symbology tab in ArcGIS Desktop 10.3
Symbology > Multiple Attributes (I am displaying species, and whether or not the tree is dead or alive) > Symbol Size > Value: DBH (diameter)

However, this menu requires that I make breaks of diameter classes to which it assigns circle symbols that do not correspond in any way to the actual diameter.

The best I can do is fiddle with the combination of breaks and ranges to get the symbol to be close to the actual diameter. This is far from perfect, and always seems to favor a certain size class (e.g. 100 cm trees are more accurately represented than 10 cm trees).

I am searching for a means of setting the circle symbol diameter based on a specific attribute parameter (in this case, tree diameter). Is this possible?


texts:
ArcGIS 10.1. When I try to export a table view with CopyRows_management, I get an invalid sql error. However, the sql statement by all accounts work fine. When I perform a GetCount_management, the results are correct. Is this a bug?

Code and error:

>>> sql = """"PARCEL_ID" IS NULL"""
>>> outTab = r"C:\Users\e1b8\Desktop\E1B8\Snapping\workspace\workspace.gdb\table"
>>> tab = r"C:\Users\e1b8\Desktop\E1B8\Snapping\Address_Tables_duplicate_addresses.gdb\_Geocode_Failures"
>>> arcpy.MakeTableView_management (tab, "tabView", sql)
<Result 'tabView'>
>>> arcpy.GetCount_management (tab).getOutput (0)
u'163083'
>>> arcpy.GetCount_management ("tabView").getOutput (0)
u'13688'
>>> arcpy.CopyRows_management ("tabView", outTab)
Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 13547, in CopyRows
raise e ExecuteError: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used.
[table] An invalid SQL statement was used.
[SELECT * FROM table] The table was not found. [table] Failed to execute (CopyRows).
>>>


texts:
I want to toggle between map view and satellite imagery in Google Earth like I can in Google Maps.

Is there an easy way to add Map View to Google Earth?


texts: Non-english Characters in cql filter in GeoServer

texts:
I have a shapefile loaded in GeoServer 2.8 and accessed by OpenLayers 2.13. This layer is in UTF-8 encoding and has some english and non-english records.
When i use english fields in cql-filter - as follows, but no differnce between WMS and WFS request- GeoServer filters data very well  and returns the result.

layer1 = new ol.layer.Image({
            title: 'Ostan',
            source: new ol.source.ImageWMS({
                url: 'http://localhost:8080/geoserver/wms',
                params: {
                    LAYERS: 'ws:layer1',
                    VERSION: '1.1.1',
                    CQL_FILTER: "NAME_ENG = 'sommy'"
                }
            })
        });

But if i try to use a field which contains non-english data in cql-filter, GeoServer returns no data. My webpage has following meta, as well:

<meta content="text/html; charset=utf-8" http-equiv="content-type" />

When i check geoserver log, i see some messy characters instead of my non-english string in cql-filter.


texts: Seeking alternative to GeoTIFF that embeds georeferencing for ArcGIS Desktop?

texts:
We have a large number of GeoTIFF images that are georeferenced. Understandably, they are extremely large, due to the TIFF format.

Ideally we would like to use jp2 images, however this image type does not allow you to embed geospatial information, like TIFF.

Is there an alternative to TIFF that allows a more compressed format, and still allows for the geospatial information to be embedded into the image (so we only have to work with a single file)?


texts: Orfeo tool box (OTB) in Ubuntu 14.04 and SAGA incompatibility (gdal versions)

texts:
I need to test segmentation of a 12000 x 12000 pixels scene.
I successfully compiled OTB in Ubuntu 14.04 and ran the process.
But when I tried to start SAGA, I got a startup error.

I found that gdal library in SAGA is version 1.10, but in OTB is 1.11.
After reinstalling Ubuntu 14.04 and SAGA with Qgis, system is running, but no segmentation in any one. SAGA exits without error message and Qgis present errors regarding image geographic extension, but is a pancromatic band.

As I see, only OTB runs very well segmentation but I need to see results and extract some polygon statistics with SAGA.
How can I install both gdal versions in parallel?


texts: Why does psycopg2 "SELECT ST_AsBinary(geom)" return ""?

texts:
I'm trying to retrieve polygons out of a postgresql database and I'm running into several problems. First, I've been able to get ST_AsText(geom) to work in conjunction with arcpy.FromWKT(), but the precision is lost as described in the PostGIS docs:

  WKT format does not maintain precision so to prevent floating truncation, use ST_AsBinary or ST_AsEWKB format for transport.

But when I try to use ST_AsBinary or ST_AsEWKB with arcpy.FromWKB I get these errors:

cur.execute("SELECT ST_AsBinary(geom) FROM master.modeled_yield_hotspots WHERE model_id = 4;")
lsTemp = cur.fetchall()
lsGeom = []
for j in range(0, len(lsTemp)):
    lsGeom.append(arcpy.FromWKB(lsTemp[j]))
arcpy.CopyFeatures_management(lsGeom, pathScratchGDB + "fcHotspots")

Traceback (most recent call last):
  File "C:\SmartFarm_Tools\python_scripts\testing.py", line 24, in <module>
    lsGeom.append(arcpy.FromWKB(lsTemp[j]))
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\__init__.py", line 1640, in FromWKB
    return gp.createObject("geometry", byte_array, None, spatial_reference)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\geoprocessing\_base.py", line 379, in createObject
    self._gp.CreateObject(*gp_fixargs(args, True)))
RuntimeError: Object: CreateObject cannot create geometry from inputs

If I print lsTemp[j] I get the following:

(<read-only buffer for 0x0000000012808190, size 82146, offset 0 at 0x0000000003092F10>,)

What is a read-only buffer? How do I get WKB out of postgis and into a feature class or shapefile?


texts: Google TMS links legal for usage in OL, GeoNode, QGIS?

texts:
I really like to use the Google tms urls in OpenLayers, QGIS (QuickMapServices), GeoNode, like https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}

But is it legal to use this links?
Will my server IP get blocked after a limit of demands?
Can i buy it somewhere or is it open for usage, unlimited usage?

I don't think that the Google Map Api is really fitting this case, because I dont use the API(no js, no styling) just the tms(xyz urls): https://developers.google.com/maps/pricing-and-plans/?hl=en

Have anybody experiences about that?


texts: Unable to create new feature for FeatureLayer with editor widget

texts:
I´m using this sample for default editor widget with the help of ArcGIS API for JavaScript. I replaced featurelayers with my own layer. It´s possible to update and delete features, but I don´t create/add new feature (console: Error: Unable to complete operation.). Can someone identify my problem ?
Can it be problem with creating of OBJECTID for new feature (I hoped that it could be created automatically) or somewhere in capabilities of service ?


texts: Field Calculator runs without error, but does not populate desired fields

texts:
I have merged numerous datasets that have the same attributes in different fields. I have TotalAcres, Total_Acres_Burned, TTL_ACRS, etc... I am trying to use field calculator to populate one field (TotalAcres) with all the acres values from the varying fields.

Using ArcMap 10.4.1 and python as parser.

Here is my code

def updateAcres ( TTL_ACRES, TotalAcres, ContrAcres, FIRE_SIZE, Field9,
Total_Acres_Burned):
  if TTL_Acres > 0:
    TotalAcres = TTL_Acres
  elif ContrAcres > 0:
    TotalAcres = ContrAcres
  elif FIRE_SIZE > 0:
    TotalAcres = FIRE_SIZE
  elif Field9 > 0:
    TotalAcres = Field9
  elif Total_Acres_Burned > 0:
    TotalAcres = Total_Acres_Burned
  else:
    return TotalAcres

Then my CodeBlock is:
TotalAcres=
  !TotalAcres!

It runs without throwing an error but it does not populate TotalAcres with the value from any of the other Acres fields.

Any ideas?


texts: Error in plotting raster extracted from an RData file in R

texts:
I am attempting to plot specific layers from a set of time-series rasters contained in RData format. I used the code below:

library(maptools)
library(maps)
library(rgdal)
library(sp)
library(spdplyr)
library(raster)
library(rgeos)

loadset <- load("MOD13Q1_MYD13Q1_EVI_97_2008_89_2009_RData.RData")
getset <- get(loadset) #successfully extracted the contents of RData file
set <- raster::stack(getset)

However when I start to plot or check the histogram using this code:

r1 <- raster(set, layer=1)

I get an error:

  Error in .local(.Object, ...) :

What does this error mean and what's the workaround?

RData files can be found HERE.


texts: Merging polygons at county level to state level?

texts:
I have a shapefile that represents counties (admin 2) only, but I need to create a state level (admin 1) polygon file (I'm using QGIS). I will need to sum up (or average) the relevant data at each district level as well.

There are far too many counties to do this manually, so does anyone know of an automated way to roll up to a higher administrative level like this?

I am new to QGIS.


texts: Update shapefile with values from Excel cell in QGIS

texts:
I have an Excel table with unique  Id values.

1
2
3
4

What I want to achieve is when editing a  point shapefile with the same Id field every time a new point is created the Id field is filled by the ids from the Excel table.
The values in the table are quite unique and have aspecific structure so I cannot use the automatic id build tools within QGIS


texts: Research study participant mapping and analysis - comparison of software

texts:
I am new to everything GIS. I have been researching different software but have found more questions than solutions.

What I want to accomplish: I'm working on a research project with participants across the USA. I have .csv files with their demographics and addresses (which I presume I can geocode rather easily). # of objects to be mapped would be in the hundreds to thousands range. Anyways I would like to create choropleth maps showing participant density in certain areas of the USA as well as heat maps or other visually pleasing displays of our participant's locations. I would then like to compare this to demographics of areas (ie race/ethnicity, age, income, etc.) which I would need to either purchase or find on the internet. This information will allow us to analyze possible reasons why certain groups choose to enroll as well as prevalence of disease in certain areas of the USA.

Budget: My company has some money to spend (if necessary - we are a non-profit) and I would like to choose the best software solution to fit our needs. I have looked into Esri/ArcGIS and MAPCITE and on the other end of things QGIS. I am trying to write a proposal comparing the possible benefits and feasibility of this kind of mapping.

Can you please give me an idea of the best option for my needs (please consider I do not have any GIS experience, but I am willing to learn).

Your guidance and wisdom are much appreciated!


texts:
I have two layers (street name and address numbers). I try to sort layers where street names covers the address numbers. Unfortunately I cannot success because address numbers will always covers street names if I change layers order also.
It is always same situation.

Why numbers (address numbers) will be prioritized over text (street names) layers?


texts: JavaScript + ArcGIS + Dropdown menu with layers

texts:
I have a script that generates a basemap with a layer that I have put on top of it. I am trying to create a dropdown menu in which the user of my HTML5 page can select the layer he wants on the basemap. I have both components: a dropdown list and a function that gives me the basemap + layer of choice.

So my question is: What is the best practice to combine these two components resulting in an interactive dropdown menu for a user to select the right layer on top of the basemap?

I have found multiple methods of getting there: http://jsfiddle.net/rn3nc/1/ and: https://stackoverflow.com/questions/1343178/change-active-li-when-clicking-a-link-jquery

Dropdown menu:

    <link rel="stylesheet" type="text/css" href="dropdown.css">

      Schaalniveau

Nederland

Gemeenten

Provincies

Wijken

Buurten

function for generating a basemap and my layer:

map.on("load", initOperationalLayer("Link to server and data"))

I want something like this dropdown menu: https://developers.arcgis.com/javascript/jssamples/dataadapterfeaturelayer.html

The code here is to difficult for me though.


texts:
I'm having a problem applying existing styles and loading new ones in ArcMap 10.1.

I go to customize > Style manager  and i see that it is empty.

After, i click on Styles tab and i select from the list the existing ESRI style and also select a new one to add.

After i click ok the Style manager appears empty as it was in the beginning , i do not know
if this is some kind of bug or either i need to enable some add-in or extension.


texts:
I'm trying to create variable with a certain function (see code below). When I run get_cloud() and print cloud from the console, it says:

NameError: name 'cloud' is not defined

When I run the for loop outside the function, it works just fine.

def get_cloud():
   for layer in iface.legendInterface().layers():
        if "cloud" in layer.name():
            cloud = layer
        else:
            pass


texts: How to style OS VectorMap district data?

texts:
Following a question about ways to load OS VectorMap data, I have imported the data into ArcGIS, but now want to style it in a OS-map style. I have searched around, and found various bits of information, but no style files to download, or instructions on using them within ArcGIS.

I've tried doing some styling manually, but the results aren't good. Is there any way to get some of the official OS styles?


texts: Converting the vertical datum of elevation data in R

texts:
I have elevation data for 699 points around Australia. I want to calculate the difference between these point elevations and the highest elevation within 25 nautical miles - in simple terms I want to get a sense of the terrain around each point.

I have been able to obtain the maximum elevation within a 25nm radius using the SRTM-derived 3 Second Digital Elevation Models Version 1.0. However, this DEM uses EGM96 as the vertical datum, whereas my point data is with reference to "mean sea level" which I take to mean AUSGeoid2020.

In order to find the difference between the two elevations I need to transform one set of data to the vertical datum of the other; it would probably be easier to transform my point data to EGM96.

I have downloaded the AUSGeoid2020 ASCII file which, I believe, will permit me to transform the point elevations to the Australian Height Datum (AHD). Now I need to convert from AHD to EGM96...

I've transformed geospatial data between horizontal datums before - in R - but I haven't been able to locate any documentation relating to transforming vertical datums.


texts: Setting Min Zoom without base map using ArcGIS API for JavaScript?

texts:
How do I set min zoom for map without defining base map(User can't zoom out beyond that label)?

Set min zoom with base map is working fine.

For example:

map = new Map("map", {
                        basemap : "streets",
                        center : [ 78.00, 21.00 ],
                        sliderPosition : "top-right",
                        maxZoom : 10,
                        minZoom: 5,
                    });

But when i am trying minZoom without base map, it's not working

basemap : ""

My main motive is use the min zoom option without base map.

Something like this in the ArcGIS API for JavaScript Sandbox.