texts: Why does map display like it is in Isometric when everything is in EPSG:4326?

texts:
I'm beginning with GeoServer and try to display a map of my country (Belgium) as an OpenLayers map (embedded into a MapPanel object). I've downloaded the municipalities limits from gadm.org and copied the file into a PostGIS database (using the tools provided in the GeoServer package). I guess it is a projection issue, but I don't get why since everything is in EPSG:4326.

Here is my code:

panel = new GeoExt.MapPanel({
   title: "MapPanel",
   renderTo: "map-id",
   layers: [
        new OpenLayers.Layer.WMS("Administrative borders",
            "/geoserver/wms",
            {layers: "bel_adm4", projection:"EPSG:4326", format: "image/png"},
            {transitionEffect:'resize', isBaseLayer:true}
   )],
   map:
      {
        projection:"EPSG:4326",
        maxExtent: new OpenLayers.Bounds(2.5412,49.5100,6.4013,51.502)
       }
 });

Here's a screenshot of the generated map:

Here's how Belgium should be:


texts: How to pass tile size of to 512 using OpenLayer 3 and GeoServer 2.14

texts:
I am trying to change the default tile size of 256 to 512 using OpenLayers 3.
Code is working fine for normal WMS request, but giving error (tiles not loading) for WMS-C request.

My sample code is:

var projExtent = ol.proj.get('EPSG:3857').getExtent();
    var startResolution = ol.extent.getWidth(projExtent) / 256;
    console.log("startResolution:", startResolution);
    var resolutions = new Array(22);
    for (var i = 0, ii = resolutions.length; i < ii; ++i) {
        resolutions[i] = startResolution / Math.pow(2, i);
    }

  var tileGrid = new ol.tilegrid.TileGrid({
        minZoom: 6,
        extent: [-13884991, 2870341, -7455066, 6338219],
        resolutions: resolutions,
        tileSize: [512, 512]
    });

   var demolayer2 = new ol.layer.Tile({
        source: new ol.source.TileWMS({

           // url: 'http://localhost:8080/geoserver/wms',
            url: 'http://localhost:8080/geoserver/gwc/service/wms',

            params: { 'LAYERS': 'topp:states', 'TILED': true, 'SRS': 'EPSG:3857'},

            serverType: 'geoserver',
            tileGrid: tileGrid

        })
    })

   var view = new ol.View({
            center: ol.proj.transform([-74.0, 40], 'EPSG:4326', 'EPSG:900913'),
            zoom: 7,
            projection: projection,
            extent: [-13884991, 2870341, -7455066, 6338219]
        });

        var map = new ol.Map({
            controls: ol.control.defaults({ attribution: false }).extend([
                new ol.control.MousePosition()

            ]),
            target: 'map',
          //  pixelRatio: 1,
            view: view
        });

        var osm = new ol.layer.Tile({
            source: new ol.source.OSM()
        });

        map.addLayer(osm);
        map.addLayer(demolayer2);

when i use url: 'http://localhost:8080/geoserver/wms' following request has been send and everything is working fine.
'http://localhost:8080/geoserver/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=topp%3Astates&TILED=true&SRS=EPSG%3A3857&WIDTH=512&HEIGHT=512&CRS=EPSG%3A3857&STYLES=&BBOX=-8875613.914302688%2C4459702.5928635085%2C-8249441.778590525%2C5085874.728575672'

But, when i use 'http://localhost:8080/geoserver/gwc/service/wms' following request has been send
'http://localhost:8080/geoserver/gwc/service/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=topp%3Astates&TILED=true&SRS=EPSG%3A3857&WIDTH=512&HEIGHT=512&CRS=EPSG%3A3857&STYLES=&BBOX=-8249441.7785905255%2C4459702.5928635085%2C-7623269.642878362%2C5085874.728575672'

and it gives following error:
'400: X,Y values for the tile index were calculated to be {18.825600646767, 39.12216711529674} which had to be rounded to {19, 39} and exceeds the threshold of 10%. Perhaps the client is using the wrong origin ?'

I have also created gridset in geoserver with 512 pixel size.


texts:
I must do a view in geoserver. I seen some example here: http://docs.geoserver.org/stable/en/user/data/database/sqlview.html. The example show the operator like for parameterizing  view. Suppose, now that in WFS request there is a value for a field (named filed1) example &viewparams=value1:a. If i use the operator like

Select filed1, filed2, field3
From  myTable
where  field1 like %value1%

the view return tuples where filed1 is a, aa, aaa, aaaa, aaaaa, etc, but i want tuples where filed1='a'. Can i use the operator =, if yes how write my query?

Select filed1, filed2, field3
From  myTable
where  field1 = .......


texts: QGIS DBmanager import file append not working

texts:
I am using QGIS and connecting to Postgresql database. I want to import shapefile data into database table via DBmanager tool.

But Append data to table selection is inactive like following image. I do not want to cerate a new database table.

How can I solve this issue?


texts: Finding or monitoring construction projects with GIS

texts:
How can we as an organization use GIS to monitor construction projects in the country, to be able to search for and identify construction taking place?

This is in order to be able to visit the various sites to ensure construction projects are registered and being built according to regulations. We are unable to manually go all over, and as a result illegal buildings are coming up that are eventually collapsing and costing lives.


texts: Changing specific dates into years in attribute table

texts:
I have a shapefile with all the buildings of the region. In the attribute table I have a field containing the date that the individual buildings have been built. I need to count how many buildings have been built in the region in certain years.

How can I change the specific dates (e.g. 15.2.2008) (the field is text form, not numeric) into year (2008)?

I use both QGIS and ArcGIS (ArcInfo), so any tips for either programme would be great.


texts: Shapefile not displayed after changing projection from EPSG:4326 to EPSG:3832 in QGIS

texts:
I downloaded the shapefile for the country "ne_10m_admin_0_countries.shp" from the natural earth web site, and opened it with QGIS and saved is with EPSG:3832.

I restarted QGIS and opened the new saved shapefile but it was not displayed and even "Zoom to Layer" did not show anything. However, I enabled on the fly projection and changed the projection to some other (I tried EPSG:4326) and it was displayed in that projection. After that check off on the fly option or change it to EPSG :3832, it was properly displayed in EPSG:3832 and "Zoom to Layer" worked properly.

Why it wasn't displayed in the first place and it was fine after I changed CRS to some other.

My environment is Ubuntu 14.04 and QGIS is 2.18.1.


texts: Including float value in loop processing of Surface Volume tool using ArcMap ModelBuilder

texts:
How do I use the Surface Volume tool in ArcMap via Model Builder, without using Python, to process a set of multiple float values in a loop,  since an iterator takes integers when I implement a for loop?
See the attached image as indication.


texts: How can I print Openlayers with base OSM and Geoserver WMS layers?

texts:
I have a map I want to print once the user open some layers.
The layers are geoserver wms and the base layer is OSM.

I know geoserver has printing module (which is installed and working because I get reply from geoserver/pdf/info.json?var=printCapabilities

What I dont understand is how I use it to print the current map with visible layers.

I've looked at GeoExt example (mapfish) but it didn't help me understand.

Update 1 :
I've generated the following URL by code :

http://77.235.53.170/geoserver/pdf/print.pdf?spec={%22units%22:%22degrees%22,%22srs%22:%22EPSG:4326%22,%22layout%22:%22A4%22,%22dpi%22:%22300%22,%22mapTitle%22:%22This%20is%20the%20map%20title%22,%22comment%22:%22This%20is%20the%20map%20comment%22,%22resourcesUrl%22:%20%22http://77.235.53.170/img%22,%22layers%22:[{%22baseURL%22:%22http://77.235.53.170/geoserver/Lehavim/wms%22,%22opacity%22:1,%22singleTile%22:true,%22type%22:%22WMS%22,%22layers%22:[%22%D7%A7%D7%95%20%D7%9E%D7%99%D7%9D%22],%22format%22:%22image/jpeg%22,%22styles%22:[%22%22]}],%22pages%22:[{%22center%22:[3875254.1134954,3680894.557955],%22scale%22:32000,%22rotation%22:0}]}

and I'm getting empty PDF (only with title)...
What could be wrong ?


texts: Which GEO data format does this data represent

texts:
I am parsing GIS data from a ESRI webservice in which I get the data specified as the following in a JSON file

"geometryType" : "esriGeometryPoint",
"spatialReference" : {
"wkid" : 102100

an example of the data provided has the following look

"geometry" : {
    "x" : 1734891.0457560429,
    "y" : 5740400.9625475956
  }

I don't usually work with GIS data so I am not sure what I am doing.
I want to store it in a SQL server DB as a GeoGrapy datatype which requires values between [-90, 90] and [-180, 180] degrees.

How Do I convert these coordinates. I do not know which format it is so I can not search for transformation equations.

could it be wgs84?

(The location is in Europe)


texts: Count points in polygons using QGIS modeler

texts:
it is the first time that I’m using the QGIS graphical modeler and I found some troubles (QGIS 2.14.1). I’m trying to count points that intersect with polygons in a vector grid, but filtering the points shapefile depending on several queries based on the attribute table. I want to do a first query and count the points, and do a second query and repeat the counting, saving both counts in the same layer. I was getting wrong results when I applied “select by expression” algorithm, even thought when I saved the selected features in an intermediate temporal shapefile before perform the counting. Using this model, I get in both counts the same results, just the one that correspond with the last selection.

Looking for in the web how to solve this problem, I found a solution in this question QGIS modeler: select by attribute wrong results , using “Extranct by attribute” instead of “select by expression” algorithm. If I’m not wrong it seems that it is a bug, so I will report that problem.

However, what I noticed is that if I perform the counting with the modeler take much more time that using the option that is in menu Vector>Analysis tools>Points in Polygons. I have read in the change log of QGIS 2.14, that PointsInPolygon processing algorithm is many magnitudes faster. It is really faster, but not when using the modeler. My question is if is there any way to use the new processing algorithm in the modeler?


texts: How do I create points in a Feature Class from a database table?

texts:
I have a destination (Feature Class):

C:/Users/XXX/Documents/ArcGIS/ABC.gdb/SamplePoints

I have a database table in the same file geodatabase that contains my source point data. I have columns with attributes in addition to the requisite coordinate data.

Using Arcpy, how do I retrieve the points from the source database table and create points in the destination Feature Class?

I tried creating a single point using this code:

import arcpy
cur = arcpy.InsertCursor(fc, ["SHAPE@XY"])
pnt = arcpy.Point()
pnt.x = 20000
pnt.y = 10000
feat = cur.newRow()
feat.shape = pnt
cur.insertRow(feat)
del cur

but, although a row is created, the coordinates are 0,0.

Current thoughts - I'll probably use a cursor to step through the records retrieved and for each row in the source, build a point and insert it into the destination. The logic seems to make sense.


texts: Moving/copying Representation Marker to Marker Symbol using ArcMap

texts:
I have a custom.style file with the marker1 Representation Marker.
Is there any way to move/copy marker1 to Marker Symbols?
(I'm using ArcGIS 10.5.1)


texts: From importing the shapefile to postgis, its geometry is added to the database postgresql. The question, how can I load the shapefile in the Visual Basic.Net using the geometry in the database ?


texts: PostGIS Shapefile Importer Projection SRID

texts:
I have a shapefile that I have imported into PostgreSQL using the PostGIS Shapefile Importer and it is currently projected in a State Plane Coordinate System. When I imported it, I selected an import SRID of 4326 (WGS84). Afterwards when I was trying to manipulate, I was getting some odd results because of the fact that my information was still seemingly in a projected coordinate system.

Does the Shapefile Importer not re-project when you assign the SRID while importing?


texts: Transforming datum after false introducing in ArcGIS

texts:
I have gps coordinates. I introduced them to arcgis without transforming them from wgs 84 to ed1950. But they are not originally ed 1950. I created a lot features with those coordinates. Then I realized the points are not at the right place. In this circumstance, what to do for not to restarting the project.

ps
My basemap has ed 1950 datum.


texts: Why can't get correct RGB img of Landsat 8 using np.dstack

texts:
I use np.dstack to stack Band4, Band3, and Band2 of Landsat8 data to get ture color image,I use code from Milos Miljkovic, but the result looks incorrect.

import numpy as np
from matplotlib import pyplot as plt
import glob
import gdalnumeric
def read_band(n):
if n in range(1,12):
    tif_list = glob.glob('*.TIF')
    band_name = 'B' + str(n) + '.TIF'
    img_id = [id for id, band_string in enumerate(tif_list) if band_name in band_string ]
    img = gdalnumeric.LoadFile(tif_list[img_id[0]])
    return img

b2 = read_band(2)
b3 = read_band(3)
b4 = read_band(4)
img_rgb = np.dstack((b4,b3,b2))
plt.imshow(img_rgb)


texts: How can I identify types of crops in agricultural raster dataset?

texts:
I have some 2014 MrSID files from the National Agricultural Imagery Program (NAIP) that I've converted to TIFF files.  I want to use them to find the area of tomato farms in a watershed.  I'm using ArcGIS 10.1.

Would an Enhanced Vegetation Index (EVI) or Normalized Difference Vegetation Index (NDVI) be useful for this?  Would I even be able to find tomato farms using orthophotos?  If there is other data out there that would be better for this task, please let me know.


texts:
I am trying to run a PCA on a set of normalized raster bands within ArcGIS. I would like to calculate the PCA on the correlation, not the covariance, matrix of the PCA, but Arc as far as I can tell only calculates PCs on the covariance matrix of the input bands. It seems the only fix to this is to feed the tool pre-normalized band values in the form of a z-score. I z-scored all of the input raster bands before feeding them into the Principle Components tool in Arc 10.1. The function is completing nicely, but the output is producing bands with mean values of anywhere between 9 and 34, even though they should be producing bands with a mean of 0, due to the nature of a PCA. What am I missing? Has anyone has experience with this issue before?

As requested, , here is a snapshot of my study site with a zscore raster overlaid: you can't see it well:

Here is a zoom-in to a New Jersey site showing PCA bands:

Finally, here is a screenshot of the first 2 PCA band output properties:


texts: How to mount an expression for "contain" in QGIS 2.18.14 Atlas/Print Composer?

texts:
EDITED: I did a refinement of the question.

Can I ask for your help? Can you download and try to help find the solution? In the link is a simplified project zip file. It has only two layers; one for roads and another for territories.- filterroads20180107.zip

I need to see in the existing Atlas / Print Composer table, filtered roads, not those that are visible on the map, but only those that belong to the current feature.

The attribute table in the 'territory' field can store more than one number. A street can belong to more than one territory.

The expression "territory" = @atlas_featurenumber ", shows the streets that have only one value, the same value of the feature, but leaves out the streets that have other numbers. e.g., show value 09, but don´t show the rows that have 09, 05, ...

1 - The expression should use a "contain" or "like", but I can not find the correct syntax. The expression "territory" like %@atlas_featurenumber%, does not work.

2- Should a "CASE" or "CASE WHEN" be used instead? How?


texts: Why does a single NAIP 2015 tile have both leaf-off and leaf-on conditions?

texts:
Image acquired at https://earthexplorer.usgs.gov/ with the following criteria:

Within tile border of leaf-off to leaf-on conditions:


texts: Running ArcPy geoprocessing tools on cursor rows?

texts:
I have a script using a Search Cursor to iteratively make copies of individual features for use as inputs in a geoprocessing tool, e.g.,

fc = <feature class>
arcpy.MakeFeatureLayer_management(fc, 'layer')
with arcpy.da.SearchCursor('layer', [<field>]) as cursor:
    for row in cursor:
           arcpy.SelectLayerByAttribute_management('layer', 'NEW_SELECTION', '{}={}'.format(<field>, row[0])
           arcpy.CopyFeatures_management('layer', 'feature')
           arcpy.<geoprocessing tool>('feature', <output>)

It runs, but is running on a large input feature class over a large raster. I thought I could speed up the script by dropping the feature copying, but can't figure out how to use a cursor row as a geoprocessing input. Is this possible?

Executing the geoprocessing tool on selected features appears to use the extent of the entire input feature class in the raster analysis, slowing the script down dramatically.


texts: Point Proximity - Extract vlaues from one point layer to another based on distance

texts:
In ArcMap I have a point layer containing two types of survey records 'Pre Survey' and 'Post Survey'. Although the Pre and Post survey locations are very similar (within 20m), they are not spatially identical. Additionally, I do not have a common unique id field linking them.

I want to be able to extract or link values from the pre and post sites.

I have looked using 'Spatial Join' which mostly works, however it seems I have to create two new point layers (a Pre and Post)? I have also looked at 'Select by Location', and 'Search Cursor' to iterate through they 'Pre Survey' records and find a 'Post Survey' record within 20m.

Is 'Spatial Join' the best way to go about this?

If 'Spatial Join' is the best way is there a way I can do this without having to create two new separate layers (i.e. can I do this from within the original point layer containing both pre and post sites)?

I have reasonable skills with scripting and python.


texts: Enable dragging on leaflet map for popup div within map

texts:
I'm sorry for my English, I just study it.
Is it possible to save map-draggable functionality on popup within leaflet map? I mean move the map while dragging(clicking and moving or tuching on mobile) across the popup?


texts: How to install OpenLayers plugin on QGIS 2.4?

texts:
I am a new user of QGIS, and I just downloaded the version 2.4 on Windows 7/32 bits. And I seem to encounter several problems:

I cannot install OpenLayers plugin
In the "Plugins" tab under "Manage and install plugins" I have only 3 windows "All", "Installed" and "Settings" and it is missing "Not installed"...and here either I open the "all" or "installed" window I have the same number of plugins, that is 23...and I cannot make any modifications like "upgrade" or "uninstall" etc...

Dos any one have an idea what is happening with my QGIS version? Or if I did something wrong when installing...I know that I just followed the easy installation for new users...


texts:
I have installed QGIS 3.2 in Windows & Ubuntu OS for the first time. 'Late comer of course'.
I cannot find the usual "standalone" QGIS Browser application installed. However, I can access QGIS Browser panel within QGIS Desktop - no problem.

Does QGIS 3 comes with 'standalone" QGIS Browser application and, if so, how do I access it?


texts:
I posted this question earlier but I like to repeat it again.

I have 50 folders. Each folder consists of 40 rasters. Each raster has a value that I need to extract. These rasters share the same x&y. Finally, I need to put them all in one data frame. To do so, I tried to run a loop to read all rasters from all folders and combine them by row. For example, I have 50 Rain rasters, so these will fall under the same column and so on. Is there any way to do that?

Here are the codes:

Reading <- function(raster) {
  df <- stack(raster)
}

All_raster <- dir("E:/50_folders/", recursive=TRUE, full.names=TRUE, pattern=".tif$")# This is to read rasters from 50 folders
All_tif <- sapply(All_raster, Reading)
ALl_tif <- as.data.frame(All_tif , xy =TRUE)

And this is the error I get
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("RasterStack", package = "raster")’ to a data.frame


texts:
"name" like 'column%'

  Output preview: 0

or

"name" = 'column_title'

  Output preview: 0

No filtering results in expression. Tried Qgis 2.12 and 2.18 both on OS X. Data is from Open Street Map - geofabrik supplier.


texts: Publishing on public site using GeoServer?

texts:
I have set up geoserver and have it running well on my LAN. I finally set up a domain and want to host my own geoserver instance publicly. I have a jetty version of geoserver running on an OS X, Mavericks Server.

What do I have to do to point my geoserver instance publicly?


texts: How to apply a gradient fill to different polygons in the same layer

texts: I am trying to make a simple map of the world in ArcMap 10.1, and I have a shapefile with all the countries in it.  I have the symbology set to unique values, so the countries are all different colors.  I want to add a gradient fill, but I can't add the gradient to its symbology without going through every single country and editing it.  If I select Properties for All Symbols, it will apply the same color to every unique value.  Any ideas?


texts: Conversion of French Lambert 1 Nord to Lat/Lon

texts:
I have used two different online coordinate converters to convert Lambert coordinates 386.1/1185.0 to Lat/Lon and receive an answer 47.42715241/-5.62147696.
This reference is 426 kms on a bearing of 242 deg from the known position of the place to which the Lambert reference is given.


texts:
I need the daily mean temperatures and precipitation values for particular days and specific localities, far from any weather gauge. I guess that, for example, this dataset kindly provided by NOAA could fit my requirements, but I would be more than thankful to anybody who may recommend me a different one with a highest resolution and easiest accessibility. I downloaded the temperature data for the specific years I need (.gz), I uncompressed them, but I was unable to do something else: I have been unable to read/load the uncompressed datasets. I have been googleing for a solution, but found nothing I was able to replicate.
Can anybody be so kind as to provide me with an idiot-proof explanation of how to handle this information in R language?


texts:
I am trying to reduce a world map raster layer (North America, Pennsylvania, etc.) I thought I could do Raster>Extraction>Clipper - I selected the raster layer and the destination of new tif file, used the select tool to select the extent area to extract. But I keep getting this error:

ERROR 1: Currently GMT export only supports 1 band datasets.

Can I only clip or extract raster layers that have 1 band only or am I doing something wrong? This raster layer has a render type of Multiband color.


texts: Cannot override onVectorElementClicked from MapEventListener

texts:
According to carto (https://carto.com/docs/carto-engine/mobile-sdk/developer-guide-to-sdk-features/), the function OnVectorElementClicked can be overwritten if it extends from MapEventListener. But I could not find it in Rider (using Xamarin (C#) Android).

Am I missing something ?


texts: Creating multiple parallel lines using ArcPy?

texts:
I'm using the following code to create multiple parallel lines

import arcpy, math

offset = 3
offset_original = offset
features=[]
arcpy.env.overwriteOutput = True

infc="Export_Output_4_SmoothLine"
talhao = "Export_Output_5"

spatialRef = arcpy.Describe(infc).spatialReference

def CopyParallel(plyP,sLength):
    part=plyP.getPart(0)
    lArray=arcpy.Array();rArray=arcpy.Array()
    for ptX in part:
        dL=plyP.measureOnLine(ptX)
        ptX0=plyP.positionAlongLine (dL-0.01).firstPoint
        ptX1=plyP.positionAlongLine (dL+0.01).firstPoint
        dX=float(ptX1.X)-float(ptX0.X)
        dY=float(ptX1.Y)-float(ptX0.Y)
        lenV=math.hypot(dX,dY)
        sX=(-dY*sLength/lenV);sY=dX*sLength/lenV
        leftP=arcpy.Point(ptX.X+sX,ptX.Y+sY)
        lArray.add(leftP)
        rightP=arcpy.Point(ptX.X-sX, ptX.Y-sY)
        rArray.add(rightP)
    array = arcpy.Array([rArray, lArray])
    section=arcpy.Polyline(array, spatialRef)
    return section

with arcpy.da.UpdateCursor(infc,["SHAPE@"]) as cursor:
    for shp in cursor:
        i=0
        while i <= 100:
            print str(i)
            if (i > 0):
                offset = offset + offset_original
            twoLines=CopyParallel(shp[0],offset)
            temp=twoLines
            features.append(twoLines)
            i=i+1

arcpy.CopyFeatures_management(features, "test")

It happens that after a while he starts to create lines in this way

How do I get these intersections or not create them in the script?


texts: How to split polygon into equal area polygons in QGIS 3?

texts:
There is a plugin in QGIS 2.18 called 'Polygon Divider' which divides an irregular polygon into smaller equal-area polygons of a specified/custom-sized but in QGIS 3.x version there is no 'Polygon Divider' plugin or tool to perform that task.

So, how to divide/split the irregular polygon into equal area/parts polygons in QGIS 3.6? Is there any Python script to split such type of polygon?

As shown in the below picture:


texts:
How we can edit the representation of a layer programmatically in DotNet ?

''Representations allow you to customize the appearance of features by storing symbol information with the feature geometry inside feature classes. This additional control can help you meet stringent cartographic specifications or simply improve the display.''


texts: Calculate Bounding Box from Center at scale

texts:
So I was trying to find my answer in Calculating bounding box by given center and scale in Android?

And below is what I've implemented (Ruby) with the size heigh/width in pixels being the same as a Letter of paper (8.5x11in @ 72dpi)

def bounding_box_from_point(center:, size:, scale_denominator:)
  dpi = 72
  inches_per_unit = 4374754
  resolution = 1 / (scale_denominator * inches_per_unit * dpi)
  half_width_deg = (size.width * resolution) / 2
  half_height_deg = (size.height * resolution) / 2

  BoundingBox.new(
    north: center.lat + half_height_deg,
    south: center.lat - half_height_deg,
    east: center.lon + half_width_deg,
    west: center.lon - half_width_deg,
    height: size.height,
    width: size.width
  )
end

Calling bounding_box_from_point(center: center, size: size, scale_denominator: scale) with:

scale = 1/18055.955520 # i.e. 0.0008861342166177423
center = Geometry::Location.new(lat: 37.806336, lon: -122.270625)
size.width = 612,
size.height = 792

It returns:

west: -122.27172131608657,
east: -122.26952868391342,
south: 37.804917238005615
north: 37.80775476199439

If you go to http://www.openstreetmap.org/export and enter those bounding box locations, you can see that the ratio does not match that of a 8.5x11in piece of paper...it's slightly too tall. What am I doing wrong here or not understanding?


texts: Is there a way to connect a GPS to ArcGIS via USB? At the moment the GPS extension only appears to support comm port connections. It would be good to know if a USB connection can be made, then perhaps my Android device can be connected and ArcGIS can use its GPS.


texts:
I've been working with GRASS GIS for about two months. I have to develop a method for the identification of target regions of small hydropower plants. Boundary conditions are a predetermined minimum and maximum flow and a minimum and maximum height of drop.
Here is a brief overview of the proposed method:

define discrete Points of River (each 50m) and give them elevation information (height)
basin determine for each Point
calculate outflow for each point and check boundary conditions
calculate: P = Q * delta_h * g * density between points

My Problem is, I don't know how to realize the first Step with Grass. I found some information about LRS (Line Reference System), v.to.point and more.

I fail to see how I put these points and correlating them height values.
What should I do so that I can make on the basis of each set point a analysis (catchment area of each point). Generally speaking, I do not know how to make an analysis that is based on set points.

My Idea was to create specific Points (e.g. each 50 meters), paste them into an attribute table (v.to.db) and then convert them to an CSV file (with v.out.ascii). Next step is to loop the r.water.outlet (example: How to do the loop when using GRASS module r.water.outlet?) for each point by using a CSV file.

However, I have my problems with it.

Do you think this is the wrong way? Your way is a completely different but sounds like a good solution. I want to understand why my approach is wrong.

Another question: How can I place equidistant vector points along a river? I can't find a good solution for that with GRASS.

The region size is about 5x5 km. Your idea is also a very good one but for my Bachelor thesis I have to do it as I describe. Maybe my method is better understood when I point my entire process here:

I also tried v.segment and v.to.points. But I can't choose an equidistant distance between points with that. So is there no other way?


texts: Field Calculator - replace string function - data value recognised as a column and is "not found"

texts:
I have the same challenge as GrahamD - How to use the replace string function in Field Calculator?

I have exactly the same challenge.

Exactly the same as GrahamD's challenge, p abbreviated for pan. Type_2 as a test column. duplicated.

Using underdark and eseglem's formulars/advice here is the copy and paste directly from the expression box:

replace("Type_2", "p", "pan")

All text and quotations in red. ok button greyed out.

  Output preview: expression invalid Parser Error: (nothing detailed) Eval Error: Column 'p' not found

I get exactly the same results should I swop Type_2 and p around.


texts: Extracting raster pixel value in ArcGIS Desktop?

texts:
I am dealing with some very large geotiff files (the whole of Canada at a 250m resolution) and is having trouble to extract the pixels value as I need to do some calculations for mapping ecosystem services (carbon,timber volume, biodiversity etc.). To do so, I need all the pixel value and we have around 4.37 billion raster pixels across Canada. I have tried Python code or convert from raster point but could not success. Then I tried to clip the raster by provenance but it only convert ~ 1/6 part of BC afterward it stopped working.

Can anyone help me to resolve this problem?


texts: How can I create a WMS-Service hyperlink to a pdf-file?

texts:
I like to create a WMS-Service with QGIS.
I want to have a link in the table of the WMS-Service where I can click on and a pdf-file opens.
I know the normal way how to hyperlink a pdf (layer properties and actions), but that doesn´t work when I use the WMS-Service.
Does anybody understand my problem and can help me?


texts:
I installed pgrouting successfully and did a little import with a 20Mb file. Works good!

Now I'm trying to import a file with the size of a country and it doesn't even start to import. The last 2 lines in the console are: "Trying to load data - Trying to parse data" and then it's done.

Why is it not working like it should? It's not even a out-of-memory bug, it's just not working!

Thanks


texts: Arcpy Insert.da.Cursor - Issues writing polygon geometry

texts:
Currently building a tool that automates the creation of some polygons that we need done quarterly.  Below is just a set script I am using that is giving me a similar issue when I try to run the script.  The current output only gives me a polygon with 4 vertices and it is not making a point for each coordinate I have listed.  This full script is much larger and complex this is just the one part I am stuck on for some reason.  This sample is one of my variations I have tried.

Code

import arcpy

workspace = r"....\Workspace\WritingGeometries"
coordinateSystem = 4326
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = coordinateSystem

fcPointTest = "TestPolygon"
fcPointPath = dsPath + "\\" + fcPointTest
arcpy.CreateFeatureclass_management(dsPath, fcPointTest, "POLYGON")

testCord3 = [arcpy.Point(-95.542806, 30.720361),
             arcpy.Point(-95.479499, 30.805696),
             arcpy.Point(-95.499301, 30.814402),
             arcpy.Point(-95.509831, 30.817513),
             arcpy.Point(-95.516297, 30.818963),
             arcpy.Point(-95.520657, 30.819754),
             arcpy.Point(-95.523791, 30.820232),
             arcpy.Point(-95.526151, 30.820543),
             arcpy.Point(-95.527991, 30.820756),
             arcpy.Point(-95.529465, 30.820909),
             arcpy.Point(-95.530674, 30.821021),
             arcpy.Point(-95.531681, 30.821107),
             arcpy.Point(-95.532535, 30.821174),
             arcpy.Point(-95.533267, 30.821227),
             arcpy.Point(-95.533901, 30.821270),
             arcpy.Point(-95.534457, 30.821305),
             arcpy.Point(-95.542806, 30.821560),
             arcpy.Point(-95.534457, 30.821305),
             arcpy.Point(-95.533901, 30.821270),
             arcpy.Point(-95.533267, 30.821227),
             arcpy.Point(-95.532535, 30.821174),
             arcpy.Point(-95.531681, 30.821107),
             arcpy.Point(-95.530674, 30.821021),
             arcpy.Point(-95.529465, 30.820909),
             arcpy.Point(-95.527991, 30.820756),
             arcpy.Point(-95.526151, 30.820543),
             arcpy.Point(-95.523791, 30.820232),
             arcpy.Point(-95.520657, 30.819754),
             arcpy.Point(-95.516297, 30.818963),
             arcpy.Point(-95.509831, 30.817513),
             arcpy.Point(-95.499301, 30.814402),
             arcpy.Point(-95.606112, 30.805696)]

#Begin Processing
array = arcpy.Array(testCord3)
polygon = arcpy.Polygon(array)
cursor = arcpy.da.InsertCursor(fcPointPath, ["SHAPE@"])
cursor.insertRow([polygon])

del cursor

Below is a screenshot of what my output is.

Solution

Once Paul helped me out and noticed an error in my coordinates I was able to get the script to run through with 0 issues.


texts: ArcGIS JavaScript API find the polygon a point belongs

texts:
I have a shapefile containing polygons and whenever a user searches an address I store latitude and longitude of the location. I want to know which polygon this address belongs and I have id assigned for each polygon in the shapefile. But I am new to ArcGIS API I couldn't figure it out. I tried this:

var queryCitiesTask = new QueryTask({
    url: "https://****/gridpointspolygon/MapServer/1"
});
var query = new Query();
query.geometry = new Point(long, lat); //
query.spatialRelationship = Query.SPATIAL_REL_INTERSECTS;
query.outFields = ["pt_v03_ind"];
query.returnGeometry = false;
queryCitiesTask.execute(query, function (result) {
    console.log("test");
});


texts: Programming Geoserver 2.0.2 to add a raster data store and layer without the UI

texts:
I have a directory of raster imagery that will be updated continually. From this imagery, I am making Image Pyramids using the Geotool's PyramidBuilder utility. I need to setup a cron job to automatically add raster datastores and layers to Geoserver without using the UI.

After looking at the REST section of the Geoserver manual I was able to add my workspace, "testWS", but trying the create an ImagePyramid datastore did not work.

Since I have access to the datastore, I expanded on the shapefile example and tried to

curl -u admin:geoserver -XPUT -H 'Content-type: text/plain' \
   -d '/opt/geoserver_data_dir/2.0.2/data/test_pyramid.pyr' \
    "http://localhost:8080/geoserver/rest/workspaces/testWS/datastores/test_pyramid.pyr/external.imagepyramid?configure=all"

Where test_pyramid.pyr is the base of my ImagePyramid at this location.

This gave me an error of "No such datastore: test_pyramid"

Is there a better way to add a new datastore and layer to Geoserver without manually adding each one via the UI? I need help in crafting the proper REST statement that will add an extisting ImagePyramid as a datastore and layer.

Is there some Java code that could do this? I looked at the Geoserver python extensions but they did not have this either.


texts: How can GeoJson polygons be colored programmatically using OpenLayers 6?

texts:
In my GeoJSON file, there is an attribute for each polygon. I am trying to update the fillcolor of the polygon base on the value of an attribute in the GeoJSON file. I have a function getColor() that returns a color, based on the input value.

There is already similar solution, but that was in 2012.

Can you help me?

  Error("${getColor(ID)}"): Unexpected template string expression. (no-template-curly-in-string)

var styles = [
  new Style({
    stroke: new Stroke({
      color: "blue",
      width: 3
    }),
    fill: new Fill({
      color: "${getColor(ID)}"
    })
  },
  {
    context:
    {

      getColor : function (feature) {
        return feature.attributes.ID > 950 ? '#800026' :
        feature.attributes.ID > 750 ? '#BD0026' :
        feature.attributes.ID > 500 ? '#E31A1C' :
        feature.attributes.ID > 250 ? '#FC4E2A' :
        feature.attributes.ID > 150 ? '#FD8D3C' :
        feature.attributes.ID > 50  ? '#FEB24C' :
        feature.attributes.ID > 10  ? '#FED976' :
        '#FFEDA0' ;
      }
    }
  })
];


texts: How does GRASS handle multi-channel images (e.g. landsat)?

texts:
I am trying to understand how to use multi-channel images (e.g. landsat or the land-cover data provided by Natural Earth). I am hoping that someone can confirm or correct my interpretation of the documentation at http://grasswiki.osgeo.org/wiki/LANDSAT and perhaps elaborate somewhat on the matter.

It appears that the appropriate module for importing multi-channel (RGB or RGBA) images is r.in.gdal. I was initially confused about why r.in.gdal split images into a raster map for each channel, but now I understand that a raster data is fundamentally a single channel of information and so the split is required if the image is to be treated as raster data. It follows then that displaying the image as a layer in a map requires re-assembling the channels, which can be accomplished with a module such as d.rgb.

I am still a little puzzled at why GRASS requires the seemingly complicated process of explicitly splitting and re-assembling the RGB channels of images. Does GRASS have some method that I haven't run across that simplifies this process to the two steps required in most other GIS software: 1) import image 2) display image as layer in map.


texts:
I'm working on trying to use arcpy to automate the creation of a report, but I'm running into a problem.

>>> arcpy.mapping.ExportReport('2013pp',"C:\reports\2013_Graffiti.rlf","C:\reports\2013_Graffiti_report.pdf")

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\mapping.py", line 506, in ExportReport
    assert isinstance(report_source, (MapDocument, TableView, Layer)), str(type(report_source))
AssertionError: <type 'str'>

I'm having trouble deciphering the traceback as it seems to reference a variety of internal python files. In trouble shooting so far, I thought the problem might be that the report_source needs to be a layer file and I was using a .shp so I converted to a .lyr but then I get the error:

Parsing error SyntaxError: invalid syntax (line 1)

In response to the suggestions below, I've tried the following with full paths, but am receiving the same error.

>>> arcpy.mapping.ExportReport('2013pp',r"C:\reports\Graffiti.rlf", r"C:\reports\2013_Graffiti_report.pdf")
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\mapping.py", line 506, in ExportReport
    assert isinstance(report_source, (MapDocument, TableView, Layer)), str(type(report_source))
AssertionError: <type 'str'>

When I try Blord's suggestion I get the following:

>>> arcpy.mapping.ExportReport(arcpy.mapping.Layer('2013pp.lyr'),"C:\reports\2013_Graffiti.rlf","C:\reports\2013_Graffiti_report.pdf")
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 389, in __init__
    super(LayerMixin, self).__init__(lyrfile)
  File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\arcobjects\_base.py", line 47, in __init__
    for arg in args))
ValueError: Object: CreateObject Layer invalid data source
>>>

The last part about the invalid data source makes me wonder about whether the shape file is acceptable...it should be given that I am able to manually create a report with it.

As an aside, is this only an arcpy mapping tool, i.e. I can't find it in ModelBuilder or in a toolbox?

I was hoping to try to trouble shoot by creating the tool in ModelBuilder, exporting it and examining the script.


texts: Removing every second row and column cells from a raster file

texts:
I have a SRTM 3 arc second topography model (.tif file).
If I remove every second row and column of cells from the grid, can it be said that the model becomes a 6 arc second?

Or does SRTM 3 arc second model needs some kind of interpolation of its grid in order to get the 6 arc second model?
If it does need, can anyone provide GDAL string example of how this can be done?


texts: Can QGIS or GRASS convert layer files to shapefiles?

texts:
Should I be downloading QGIS or GRASS?

I want the program with the most features. Most importantly, I have a lot of .lyr files. Can either of these softwares convert to shp and if so, how can I do it?

I've tried exporting through grass, it runs but no image appears. It says something about "try setting parameters".


texts: Configuring legend for ArcGIS API for Silverlight?

texts:
I can't get the simple Legend sample working with this florida mapservice and the ESRI Silverlight SDK 2.1.

The names appear, but the tree doesn't get populated.

Anyone see why this code shouldn't work?

Here's my modified xaml:

<esri:Map x:Name="MyMap" Extent="-15000000,2000000,-7000000,8000000">
            <esri:ArcGISTiledMapServiceLayer ID="Street Map"
                    Url="http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"/>
            <esri:ArcGISDynamicMapServiceLayer ID="Florida" Opacity="0.6" Initialized="ArcGISDynamicMapServiceLayer_Initialized"
                                               InitializationFailed="ArcGISDynamicMapServiceLayer_InitializationFailed"
                    Url="http://map.floridadisaster.org/GIS/rest/services/Facilities/Critical_Facilities/MapServer"/>
        </esri:Map>

        <Border Background="#77919191" BorderThickness="1" CornerRadius="5"
            HorizontalAlignment="Right"  VerticalAlignment="Top"
            Margin="20" Padding="5" BorderBrush="Black" >
            <esri:Legend Map="{Binding ElementName=MyMap}"
                         LayerIDs="Florida"
                         LayerItemsMode="Tree"
                         LayoutUpdated="Legend_LayoutUpdated"
                   BindingValidationError="Legend_BindingValidationError"       />
        </Border>

Running the sample with fiddler, I see a request sent to here:
http://www.arcgis.com/sharing/tools/legend?soapUrl=http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Specialty/ESRI_StateCityHighway_USA/MapServer&f=json&returnbytes=true&_ts=634254399369803267

However, no similar request is made when running my modified version.


texts:
I have a line which I intersected with a polygon layer in order to split the line along the polygon borders and add the polygons' attributes to the lines. Further processing of this data happens in Excel - but before this I need to add an ascending ID along the line for easier orientation.
I am able to add the current $ID of the lines, but it is not sorted/ascending.

How can I add an ascending ID which follows the line direction?
Using QGIS 3.10


texts:
I have to update a rank column from 1-5  based on the category field.
For example if my category is COMGOV assign rank 1
If category is HOTPRE OR RCNCLB or HOTRES assign rank 2
If category is COMBAN or HLTSPA or so on assign rank I need to run this using arcpy script...

import arcpy
Tl= "D:/Shreya/DL.shp"
fields = ["FTR_CRY", "RANK"]
with arcpy.da.UpdateCursor(Tl, fields) as cursor:
   for row in cursor:
      if row[0] =="COMGOV":
         row[1] = 1
      cursor.updateRow(row)


texts: Why does Delimited Text Plugin display red sample columns?

texts:
When importing a CSV file (with associated CSVT) using the Add Delimited Text Layer feature, sometimes the preview window displays a column of sample text in a red font.

Does anyone know the reason for this? I thought it might be a data type mismatch with the CSVT file but everything seems fine with it.


texts: How to "Merge Selected Features" with python?

texts:
I have been looking for the answer to this question for some time now, but can't seem to find it! My problem is straight forward:
I have one layer
It contains lots of features (only polygons), say countrys
I want to merge some features together.
Using python.

I know that I can select features with the mouse and hit the "Merge selected features"-button, but that is to slow.
I can use python to select the features that i want, but how can i merge the selected features from python? (the features don't have a "common field" or something like that, what features are to be merged will be decided by in-data)

Help is greatly appreciated!


texts: Is it possible to place a point at a coordinate in QGIS?

texts:
Is it possible to place a point at a coordinate in qgis?

I have a map, using irish grid reference and wish to place a point at an exact 6 digit coordinate. Is it possible to do it.

Previous to this i have placed the point roughly and then moved it whilst zoomed to the exact location, surely there is an easier way.


texts: Writing masked GeoTIFF in rasterio and visualizing in QGIS?

texts:
I have an interpolated raster, which is a masked array, that shows where the data is valid. I want to save it as masked TIFF with rasterio.

To test, here's the JSON file, and the code:

import rasterio, rasterio.io, rasterio.transform
import numpy as np
import json

with open('data.json') as f:
    d, m = json.load(f)

zi = np.ma.masked_array(d, m)
res = 100000
minx = 2282862
miny = 5309393
dtype = 'uint8'
tr = rasterio.transform.from_origin(minx - res / 2, miny - res / 2, res, -res)

raster_params = {
    'driver': 'GTiff',
    'width': zi.shape[1],
    'height': zi.shape[0],
    'count': 1,
    'dtype': dtype,
    'compress': 'LZMA',
    'crs': {'init': 'epsg:3857'},
    'transform': tr
}

with rasterio.Env(GDAL_TIFF_INTERNAL_MASK=True):
    with rasterio.open('raster.tif', 'w', **raster_params) as dataset:
        dataset.write(np.array([zi.astype(dtype)]))
        dataset.write_mask(zi.mask)

As a result, there's one band and it has zeroes instead of nan and no mask.

I tried using two layers and writing this way:

...
    'count': 2,
...
        dataset.write(np.array([zi.astype(dtype), zi.mask.astype(dtype)]))

It creates a band that I can choose in QGIS. But it just hides the layer completely. When Band 2 has '1' for data and '0' for no-data, it looks like this:

If I do the other way around ('1' for no-data, '0' for data), the result is the same.

Am I writing it correctly to GeoTIFF?

How to visualize it properly in QGIS?


texts: Convert raster to polygons, crop, reproject

texts:
I wish to map a high-resolution grid of population counts of Korea and extract from it the Seoul Capital Area (SCA) at the same level of resolution.
The population data comes in two forms: raster (population_kor_2018-10-01_geotiff.zip) or CSV centroid data (population_kor_2018-10-01.csv.zip). The SCA administrative boundaries come from GADM (R (sf) file).

Am I doing this right? in the right order? Is it not correct to reproject the raster before converting to polygons since the variable (population) is areal thus should be recalculated? I have also seen it recommended to reproject polygons rather than the raster.

library(raster)
library(rgdal)
library(sf)

#Load raster of population
kor_dens <- raster("population_kor_2018-10-01.tif")
kor_dens@crs #CRS is WGS84

#Load SCA admin boundaries
sca_gadm <- readRDS("gadm36_KOR_1_sf.rds")
st_crs(sca_gadm)$epsg #CRS is WGS84

#Crop raster down to SCA first to save on computing (rectangular outcome)
kor_dens <- crop(kor_dens, extent(sca_gadm))

#Reproject raster to planar CRS, which adjusts population values
kor_dens <- projectRaster(pop_rast, crs=CRS('+init=EPSG:5179'))

#Reproject SCA boundaries to planar CRS
sca_gadm <- st_transform(sca_gadm, 5179)

#Mask raster, which replaces values outside the area of interest with NA
kor_dens <- mask(kor_dens, as(sca_gadm$geometry, 'Spatial'))

#Convert raster to polygons, where cells with NA are not converted
kor_dens <- rasterToPolygons(kor_dens)

The resulting polygon object is massive at 9.3 Gb, so I'm still not sure this is the way to go for statistical analysis.

EDIT: I've replaced the code entirely, inspired by @BrunoConteLeite's answer.


texts:
I have a problem with integer types when importing a shapefile created in OGR into PostGIS.

I started off with a field created as ogr.OFTInteger64 in OGR/pyqgis as follows:

field_defn = ogr.FieldDefn( fld, ogr.OFTInteger64 )
lyr.CreateField ( field_defn )

The field type shows as type: qlonglong, and type name: Integer64 in QGIS.

But after importing the shapefile into PostGIS with shp2pgsql,

shp2pgsql -Id -s 4326 my_lyr.shp my_lyr | psql

the field ended up being a double precision field in \d+, and I had to convert values in this field back to BIGINT before use.

According to PostgreSQL documentation, both bigint and double precision are types 8 byte (64 bit). How can I avoid the confusion and get the intended BIGINT field in PostGIS?


texts: Export image collection Google Earth Engine to Google Cloud Bucket Python API

texts:
I have an image collection within GEE's Python API that I want to clip and write to my google cloud bucket. I am working on a Google Cloud VM and have mounted my bucket "gee-bucket22" with gcsfuse. I have successfully written images from my ee.batch export in the past, but for an unknown reason, my task is not executing. The task does not report an error and does not successfully write images. The status of the task is 'state': 'READY' when calling task.status(). Is there an additional step to execute tasks I'm missing?

import ee
ee.Initialize()
AOI = ee.Geometry.Rectangle([115.647,-33.669,115.690,-33.685])

imageCollection = ee.ImageCollection('MODIS/006/MOD13Q1').filterDate('2017-01- 01', '2018-10-01').filterBounds(AOI).select('EVI')
collectionList = imageCollection.toList(imageCollection.size())

out_image = ee.Image(collectionList.get(1))#export first image
task = ee.batch.Export.image.toCloudStorage(
       image=out_image,
       bucket='gee-bucket22',
       fileNamePrefix='exampleExport')
task.start()
print(task.status())

I'm sure my bucket is mounted because after running

gcsfuse gee-bucket22 ~/mnt/bucket-location

I receive a response "bucket successfully mounted".


texts:
I have a year worth of daily snow cover data in tif files.
I want to merge them all into one (stack) and select a polygon from which I want to extract snow cover extent as in sqare miles and as % of total area.

I want to create time series of all the daily files for a selected square location area.


texts: Creating Ellipse in WGS84, but with metric parameters

texts:
I have the function

--Ellipse(x,y,rx,ry,rotation,#of segments in 1/4 of ellipse)
SELECT ST_Translate( ST_Rotate( ST_Scale( ST_Buffer(ST_Point(0,0), 0.5, $6), $3, $4), $5), $1, $2)

I have points in WGS84 (x,y), radiuses in meters (rx, ry), rotation angle in degrees. And I need a new WGS84 geometry as a result.

As you can see I can't call this function without conversions between WGS84 and some metric SRID. But my points aren't in a definite region. So a SRID which is optimal for one point, throws an error for another point during conversion.

Is there a common way to get an optimal metric projection (SRID) for a WGS84 point? Or at least to use some rough, but world-wide universal metric projection? Or maybe there's a trick for this case?


texts: Source for countries/country codes for the entire world?

texts: Is there any Google API which would return all the countries in the world with the country code?


texts: SQL expression for finding strings including letters in ArcGIS Desktop?

texts:
I'm working on ArcGIS 9.3.1. I have a shapefile with a string field containing point names like these:

10
12
12/13
t/9
6577y
1234
12t9
124y67
u9

and so on (and yes, I know it looks like a colossal mess. It's not my design). Using  WHERE MyColumn LIKE '%Search String%' I can find (using select by attribute) the ones containing a specific letter or a backslash (/).

How can I select all values containing a letter (in this example: t/9, 6577y, 12t9, 124y67, u9), not caring what letter it is or it's position? (is there a way to specify a position?)


texts: applying a gradient to terrain curves in QGIS: problem

texts:
I'm trying to apply a gradient to terrain curves generated from a geotiff. I then choose a graduated style with a gradient color palette.

everything seemed fine until I noticed that from a certain altitude, color is wrong. look at the top of peaks, suddenly light-blue (= low altitude) instead of dark blue (= high altitude)

geotiff looks fine:

Anybody knows the problem or has a solution?


texts: Required version of proj4.js with OpenLayers 4.6.4

texts: I am currently using the version 2.2.2 of proj4.js with OpenLayers 4.6.4 and I plan to move to the version 2.4.3 of proj4.js because it resolves a defect that I have. However, I wanted to know if OpenLayers 4.6.4 works with proj4.js 2.4.3.


texts: Using NetCDF for point time-series observations?

texts:
I have a large growing dataset of water and climate observations for several years in regular hourly time-step.
The number of locations is cca. 1000 locations with latitude, longitude coordinates.
At each location up to 5 parameters (water level, discharge, rainfall, snow depth, temperature) are measured at regular hourly time-step (gaps in time-series are labeled with a special "nodatavalue" code).

Currently I am using a Microsoft SQL Server database with a following main tables:

stations (id, name, latitude, longitude, elevation)
parameters(id, name, units, scalefactor, nodatavalue)
stations_parameters (station_id, parameter_id, start_time, end_time)
waterlevel (station_id, time_utc, value)
discharge (station_id, time_utc, value)
temperature (station_id, time_utc, value)
rainfall (station_id, time_utc, value)
snow (station_id, time_utc, value)

The size of my database is 2 GB. For various reasons mainly due to database size limitation on my webhost and interoperability issues I need to move away from MS SQL Server. I would like to use the NetCDF format because I read that it is suitable for multi-dimensional data, and that it has built-in support for space and time queries. Specifically I need to run the following types of queries very fast:

Query 1: For a time-range give me the average, max, min or sum of a parameter for all stations. For example create a map of maximum temperature between 2013-07-25 and 2013-07-26.

Query 2: For a time-range, point location and parameter, give me the time-series. For example hourly time-series of temperature at station 777 from 2013-01-01 to 2013-07-25.

My questions:

Is it possible to use NetCDF for this type of multi-dimensional
data? (note that the space is not a regular grid but rather sparse
point-locations)
Is there any free software tool that has built-in support for visualizing a point time-series NetCDF data file? I have used
Panoply and IDV but only with regular NetCDF grids.
If yes - What would be the structure of my NetCDF data file (dimensions, variables and their order)?


texts: How to add multiple map marker while reading data from a single csv file

texts:
I am trying to add a different map marker for every lat/long point added through a csv file using CSV Layer,is it possible to loop through each of the point for which I am going to add map marker and add a different map marker for those points.

I am having a CSV file that contains lat/long and other information such as device_type for various electronic devices such as below:

latitude,longitude,device_type

28.7041,77.1025,device_1

26.8541,76.8965,device_1

27.7844,79.7974,device_2

26.7041,76.1025,device_2

On the basis of device type,I want to create a different map marker and plot it on map,like for the device_1,I want to create one type of map marker and for the device_2,I want to have a different map marker.Is it possible to do so?

    addCSVLayer = function (){
        // Importing csv file
            var url =
            "http://localhost/Geocoding_Report.csv";
            var template = {
            title: "Device Info",
            content: "Device type {device_type}"
            };

            csvLayer = new CSVLayer({
                url: url,
                popupTemplate: template,
            });

            csvLayer.renderer = {
                type: "simple", // autocasts as new SimpleRenderer()
                symbol: {
                  type: "picture-marker",  // autocasts as new PictureMarkerSymbol()
                  url: "http://localhost/fibrenode.png",
                  width: "25px",
                  height: "25px"
                }
              };
            map.add(csvLayer);
        }
    }


texts: Export ArcGIS Online feature service with attachments

texts:
How can one export an attachment enabled feature service from ArcGIS Online and maintain it's integrity, meaning: get the geometry, table values, attachments, and the relationship class which connects them?

The route of "Open in desktop" > select layer > Data > Export only saves the geometry and table values. Ditto for Feature Class to Feature Class tool.


texts: Seeking ArcGIS Pro tool to batch change data sources?

texts:
ArcCatalog has a useful Set Data Source tool which allows for rapidly updating multiple data sources in a map document that have changed.

Is there a similar functionality in ArcGIS Pro for project files, or within a project?

I can presumably also script a solution like Repairing data sources in multiple MXDs using ArcPy?, but not all our team members have the same skill level.


texts: Cannot convert dms to dd in Excel 2007 using VBA code

texts:
I am trying to convert the data in an Excel sheet from dms (degree-minute-second) to dd (decimal degree) format using the vba script available here.

Here is sample data:

53° 0'  39.31"  6° 44'  50.14"
52° 59'  51.43" 6° 46'   0.02"

In the Excel sheet, this code is working:

=Convert_Decimal("10° 27' 36""")

but this is not working:

=Convert_Decimal(B3).

The error is #value!

I want rows in excel in dd format. Alternatively I want to convert this Excel file into a database but the format must be dd.

Any suggestions anyone?


texts:
I've been using ColorTool (based on supreme ColorBrewer) for most of my mapping needs in ArcMap. After migrating to ArcGIS 10 I discovered that it will no longer be a case.

Are there any other options to include ColorBrewer palettes in ArcGIS 10?


texts:
I am trying to resample a landsat NDSI image to a MODIS NDSI image. When I run the code below and use the inspector to click on a few random points in the image, it provides the normalized difference value at the new resolution, but then in parenthesis it provides a percentage. This not at all the points, but it looks like it is at most. I can't seem to find any info on this in the GEE developers guide. What does the percentage value correspond to? And how should I interpret it?
// Get information about the MODIS projection.
var Vis = {
  palette: ['000088', '0000FF', '8888FF', 'FFFFFF'],
};

var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD10A1')
    .filter(ee.Filter.date('2018-10-28'))
    .first())
    .select('NDSI_Snow_Cover');

var modisProjection = modis.projection();
print('MODIS projection:', modisProjection);

// Load and display forest cover data at 30 meters resolution.
var l8NDSI = ee.Image('LANDSAT/LC08/C01/T1_RT_TOA/LC08_200016_20181028')

// Load another image to replace the cloudy pixels.
var cloudMask = ee.Algorithms.Landsat.simpleCloudScore(l8NDSI).select('cloud').lte(20)
var l8NDSInoClouds = l8NDSI.updateMask(cloudMask)
var l8NDSIsnow = l8NDSInoClouds.normalizedDifference(['B3', 'B6'])

//////
var snow = l8NDSIsnow.select('nd')

Map.centerObject(snow,10)
Map.addLayer(snow, Vis, 'snow cover 30 m');
//Map.addLayer(l8NDSInoClouds, {bands: ['B4', 'B3', 'B2'], max: 0.4}, 'masked');

// Get the snow cover data at MODIS scale and projection.
var snowMean = snow
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

// Display the aggregated, reprojected forest cover data.
Map.addLayer(snowMean, Vis, 'snow cover at MODIS scale');


texts:
I am in process of migating a web ADF 931 application to WEBADF 10. The application uses some samples from ESRI forum/samples. I am curious if there are versions of these files available within WebADF 10 environment. it seems like these are designed to work with 931 web controls. Appreciate any feedback. Thanks Jay

Referenced libraries:

<%@ Register Assembly="gridResults, Version=1.3.0.0, Culture=neutral, PublicKeyToken=4117b0729e6ae73a" Namespace="esri_samples" TagPrefix="esri_samples" %>

<%@ Register Assembly="PrintTask_VBNet, Version=2.0.0.0, Culture=neutral, PublicKeyToken=2a7e7966c16048d7" Namespace="PrintTask_VBNet" TagPrefix="printTaskVB" %>

<%@ Register Assembly="WorldViewSolutions.ArcGISServer.WebControls.FeatureSearch"
Namespace="WorldViewSolutions.ArcGISServer.WebControls.FeatureSearch" TagPrefix="wvs2" %>

These are the errors I am getting:

Unable to create type 'esri_samples:gridResults'. Could not load file or assembly 'ESRI.ArcGIS.ADF.Web.UI.WebControls, Version=9.3.1.3000, Culture=neutral, PublicKeyToken=8fc3cc631e44ad86' or one of its dependencies. The system cannot find the file specified.

Unable to create type 'wvs2:FeatureSearchTask'. Could not load file or assembly 'ESRI.ArcGIS.ADF.Web.UI.WebControls, Version=9.3.1.3000, Culture=neutral, PublicKeyToken=8fc3cc631e44ad86' or one of its dependencies. The system cannot find the file specified.

Unable to create type 'oat:OpenAttributeTableTask'. Could not load file or assembly 'ESRI.ArcGIS.ADF.Web.UI.WebControls, Version=9.3.1.3000, Culture=neutral, PublicKeyToken=8fc3cc631e44ad86' or one of its dependencies. The system cannot find the file specified.


texts: What files needed to merge two shapefiles using QGIS?

texts:
I am trying to build a choropleth map for an assignment. This requires that I merge the shape files in the following compressed folders into a single shapefile.

folder 1
folder 2

I have just installed QGIS to get this done. My challenge is that I like to know the contents in these folders that are needed and those I should discard.
I am new to GIS.


texts: Export pixel neighborhood in Google Earth Engine

texts: I have some points and know how to extract the Landsat data at these points and build and export a training dataset to feed an external classifier. But how can I efficiently and automatically build a small patch of 2-D pixels (say 3x3 neighborhood) around each sample point to build and export a similar training dataset to feed an external convolutional classifier?


texts: Consistent color scheme across layers of different types

texts:
I'm trying to figure out a way to get a consistent color scheme across layers of different types.

I have a layer of lines (electrical cables) and a layer of points (electrical nodes). The elements in those layers have an integer attribute that corresponds to the ID of their electrical power source.

Simply put, I'd like the cables and node that share the same power source to have the same (random) color.

So I made a categorized style with random colors on one of the layer but I cannot copy it to the other layer since I cannot copy a style from a layer of points to a layer of lines. And if I make another categorized style on the other layer, the random colors will be different.

Also, I can not copy each color by hand since I have about 200 different possible values of the attribute.

Is there a reasonably simple way to achieve what I want to do ?


texts:
I have code like this from Xander Bakker

def main():
    import arcpy
    fc = r'Drive:\Path\To\Source\FC'
    fld_seq = 'SequenceLbl'

    # Add field if it does not exist already
    AddField(fc, fld_seq, "TEXT", 12)

    # start update cursor to update the values in the new field
    flds = (fld_seq)
    cnt = 0
    with arcpy.da.UpdateCursor(fc, flds) as curs:
        for row in curs:
            cnt += 1
            km, hect = divmod(cnt, 10)
            hect_txt = "%03d" % (hect * 100,)
            lbl = "{}+{}".format(km, hect_txt)
            print cnt, hect_txt
            curs.updateRow((lbl, ))

def AddField(tbl, fld_name, fld_type, fld_length):
    if len(arcpy.ListFields(tbl, fld_name)) == 0:
        arcpy.AddField_management(tbl, fld_name, fld_type, None, None, fld_length)

if __name__ == '__main__':
    main()

the result will show
0+100
0+200
0+300
..
0+900
1+000
1+100
and so on..

I want to start from 0+000 instead of 0+100,
what code do I need to add or modify from that code?


texts:
I have several jpeg2000 tiles non georeferenced and covering a certain area. I know in qgis you can mosaic image using the merge tool but I have done it only using georeferenced images, which is quite straightforward.
I would like to know if there is a method with qgis to mosaic raster images non goreferenced. So basically just stitching tiles above, below, left and right.


texts:
I have a Leaflet map that I show streets and also there is an unordered list of street names in the same page with map. What I want to do is when mouse cursor hovers on a street name in the list, the poly-line which belongs to that street should be bold.

As far as I know, there is no method listed for Leaflet like that by default but there are options for poly-lines for styling and mine could be shown as below:

var polyline = L.polyline(latlngs, {weight: 6}).addTo(map);

How can I apply the bold styling to my poly-lines on a hover event?


texts: Identifying features using ArcGIS Engine?

texts:
I am new to arcgis.

I am building a custom desktop GIS application using arcgis engine 9.3.1 .

I want to make a custom identify tool and add it to toolbar control.

Any example or code will be helpfull.


texts:
I have a D8 flow direction raster made from with TauDEM's D8FlowDir method and a point vector with an attribute 'id' for each.  For some reason when I use these two as input in the GageWatershed method I get a single large watershed instead of a watershed for each point.  I checked several pages to make sure my point vector has the right attribute name.  Could someone please advise.

I have the D8 raster, point vector, and bad result watershed raster located at the following:

ftp://ftpext.usgs.gov/pub/er/va/reston/gpetrochenkov/


texts: How to replicate an unversioned SDE geodatabase to a file geodatabase?

texts:
I need to implement a process for creating a file geodatabase copy (a read-only snapshot, essentially) of a production SDE geodatabase without using versioning, adding global ID columns, etc. I cannot change any aspect of the SDE geodatabase.

All of the documentation I have read on geodatabase replication says that the parent replica data must be versioned. Indeed, attempting to create a simple check-out replica with either the Create Replica wizard in ArcMap or the Create Replica geoprocessing tool fails with the following error message:

  Can't create empty replica. Read only data or unversioned data cannot
  be replicated. Data versioned with the option to move edits to base
  cannot be replicated. GlobalIDs are required for two way and one way
  replica data. Creating one way archiving replica requires archiving to
  be enabled on the parent.

Is there a way around this?

I need to be able to run the process once a month so a scripted or programmatic solution would be preferred (arcpy, ArcObjects., etc.). The amount of data to be synchronized is large, probably 500 datasets/tables/feature classes.

One solution that comes to mind is exporting and importing ESRI XML Workspace documents, but this seems heavy-handed and error-prone. Another would be a more fine-grained approach of listing and copying individual datasets with arcpy or ArcObjects but if there is an easier, existing solution (i.e. geodatabase replication, if it would let me), that would be preferable.


texts: import parameters for r.watershed - (watershed basin analysis program) tool from grass using Qgis 2.2

texts:
I want to run this tool from Qgis 2.2 but i don't know how to use the PARAMETERS to run this tool.

The red frame are the PARAMETERS I don't know.

Can someone to tell me what those PARAMETERS are and what they mean?


texts: Field Calculator expression to replace string only if it occurs

texts:
I have two fields (STREET1 & AKA_STREET_TEST1).

I'm having a problem creating a simple IF THEN Python statement in the Field Calculator box.  I've been searching the threads for a similar example but to no avail.  Basically the formula would work something like this.

  If STREET1 contains "AVENUE", calculate AKA_STREET_TEST1 to "AVE",
  ELSE leave it blank

(example below)

Below is my feeble attempt at creating the if-then statement.


texts:
I am creating an action in QGIS, which opens a link from an attribute.

Is it possible to open such a link in QField (Android App)? Or is there any other way to open links from the attribute table in QField?


texts: How to use ATOM FEED Data in QGIS 3.6.1?

texts:
I would like to know how I can work with ATOM FEED Data in QGIS. In ealier QGIS versions this  worked with the Plugin "INSPIRE ATOM Client", but this Plugin doesn't work in the current version.

  An Atom feed is a web feed that uses the Atom Syndication Format and the Atom Publishing Protocol. Atom comprises an XML-based format that describes an Atom feed and the items of information in it, and a protocol for publishing and editing Atom feeds. -https://www.ibm.com/support/knowledgecenter/en/SSGMCP_5.5.0/fundamentals/web/dfhtl_atom_whatis.html


texts: GroundFilter.bat creating a .lda file and a 1 KB laz file


texts: Creating Polygons using Google Earth Pro?

texts:
I am brand new to Google Earth Pro, and in walking through the tutorials, my polygons don't work as they're supposed to.  I follow all of the instructions from the Google Earth lessons, but nothing shows up on my earth.  The steps I follow:

Fly to address/location;
Click the +Polygon button;
Enter name, choose colors & density, make altitude Relative to Ground at 5 meters, click OK.
Get sad that nothing happens.

I've tried with several different locations, but I get the same result every time.

Am I not filling in something that is required?


texts:
I have been trying to convert a large OSM file (5G) to the GeoJSON format. All the tools I have been using, use extensive memory and throw an out of memory exception.

What is the best way to convert a large OSM file to the GeoJSON format? Do I have to split the OSM file?


texts: Casting idataset to ifeaturelayer in ArcObjects C#?

texts:
I am try to casting idataset to ifeaturelayer i did not get any idea

here my code

            ESRI.ArcGIS.ArcMapUI.IMxDocument mxd = ArcMap.Application.Document as
                                                   ESRI.ArcGIS.ArcMapUI.IMxDocument;
            ESRI.ArcGIS.Carto.IMap map = mxd.FocusMap;
            ESRI.ArcGIS.Geodatabase.IWorkspaceFactory wsf = new
                        ESRI.ArcGIS.DataSourcesGDB.FileGDBWorkspaceFactory();
            ESRI.ArcGIS.Geodatabase.IWorkspace ws = wsf.OpenFromFile(@"D:\ARC OBJECTS\TemplateData.gdb", ArcMap.Application.hWnd);
            ESRI.ArcGIS.Geodatabase.IFeatureWorkspace fws = ws as ESRI.ArcGIS.Geodatabase.IFeatureWorkspace;
            ESRI.ArcGIS.Geodatabase.IFeatureDataset fd = fws.OpenFeatureDataset("World");


texts:
I have created a uniform (nominally square) gird of polygons using the grid tool.
I have a road network layer that consists of lines. What I would like to do is iterate over the grid polygons and clip the road network within each square.
I can certainly clip the road network out using a single polygon but I would like to automate the process :-) (160+ polygons in the grid!).

Happy to use QGIS 2.14 or 2.18!


texts: Extract many or clustered cells attributes around points using SAGA GIS?

texts:
My question is similar to Extracting raster values at points using Open Source GIS?, but perhaps more complex.

How to extract the attributes of clustered cells around a point/points?

For example, I have many vegetation plots in 100m x 100m and have to examine the spectra at the points over SPOT Data (satellite image) which has 10 m resolution.


texts:
I am running: gdal_contour -a elev hazard_interp4.asc haz_ak_contour_4_10 -i 10

the hazard_interp4.asc file has values in band 1 (singel band gray).

With the above command, I get a shapefile that has no values.

If I burn a color schema into the raster with gdaldem color-relief, then I can successfully create a contour shapefile. However I would prefer to not have to do this step.

Update:
the output of hazard_interp4.asc is:

Driver: AAIGrid/Arc/Info ASCII Grid
Files: hazard_interp4.asc
       hazard_interp4.asc.aux.xml
Size is 300, 300
Coordinate System is `'
Origin = (-19003400.000000000000000,11484799.899999998509884)
Pixel Size = (10698.000000000000000,-15086.032999999999447)
Corner Coordinates:
Upper Left  (-19003400.000,11484799.900)
Lower Left  (-19003400.000, 6958990.000)
Upper Right (-15794000.000,11484799.900)
Lower Right (-15794000.000, 6958990.000)
Center      (-17398700.000, 9221894.950)
Band 1 Block=300x1 Type=Float32, ColorInterp=Undefined
  Min=0.005 Max=0.595
  Minimum=0.005, Maximum=0.595, Mean=0.141, StdDev=0.141
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=0.59498530626297
    STATISTICS_MEAN=0.14081028623394
    STATISTICS_MINIMUM=0.0050000008195639
    STATISTICS_STDDEV=0.141377304974


texts: How to display a line with fading transparency in qgis?

texts:
I have a straight line with two vertices, which I'd like to style. I want the style to be a fading line from start to end, so the transparency is gradually increasing from start to the end of the line.
I thought I could make it with geometry generator, using some function to locate point along the line, then splitting the line into smaller parts and applying something like color_hsva() with scale_exp() to map the alpha to a distance. Unfortunately, I couldn't find any function to get a line substring or any analogue to the st_line_locate_point() postgis function.
Of course I can do this by splitting line in the data, adding a segment number to each part from startpoint, and then style them separately, but it's a bad solution.
Any ideas how to do this by qgis styling and without messing with the data?


texts: ESRI Print Task (export web map) returns 500

texts:
I tried to export web map(Esri) with print task service. it works in localhost and in staging.xxx.com(http) too.But in preview.xxx.com(https) , it throws the error 500: Error performing execute operation.

{"error":{"code":500,"message":"Error performing execute operation","details":[]}}

I can't figure out what's happening , staging uses http and preview uses https , that's the only difference.

Here's my code

        var printTask = new esriLoader.PrintTask(PrintTaskServiceUrl);
        esriLoader.Config.defaults.io.proxyUrl = proxyUrl;
        esriLoader.Config.defaults.io.alwaysUseProxy = true;
        var template = new esriLoader.PrintTemplate();
        template.exportOptions = {
            width: 1015,
            height: 633,
            dpi: 96 // if 200 ,map image will get good quality
        };
        template.format = "JPG";
        template.layout = "MAP_ONLY",
        template.preserveScale = false;
        template.showLabels = true;
        template.showAttribution = false;

        template.layoutOptions = {
            "legendLayers": [], // empty array means no legend
            "scalebarUnit": "Miles",
            "copyrightText": "
xxx
",
        }
        var params = new esriLoader.PrintParameters();

        params.map = map;
        params.template = template;

        printTask.execute(params, success, error);


texts: How to identify and locate maximum values in each individual row of a raster?

texts:
I am looking to identify and locate maximum values in each individual row of a raster. I have tried many different methods including edge detection (Sobel, Canny,   Prewitt, etc.), focal statistics and generating contours (ArcGIS). However, the magnitude of the maximum value vary considerably from row to row - so the lower value maximums are difficult to (visually) identify using the methods mentioned above (see sample image below).

Ideally I would like to start at the top row of the raster and identify the location of the maximum value in each row from top to bottom.


texts: Leaflet: how can I update content of popup by jQuery

texts:
in my popup is a spinner with class=spinner and an empty container for addon informations.

I want to hide after successful ajax.
I want to add some html in a container.

I tried this:

~~~

    map.on('popupopen', function(e) {
        var marker = e.popup._source;
        $.getJSON(url,function(data) {
        $(e.popup).find('.spinner').hide();
        $(e.popup).find('container').html('same stuff')

});
~~~

It doesn't work. Any ideas?


texts: Empty geometry resulting from pgrouting function

texts:
Let me start by saying I'm assuming I'm missing something very silly but I can't see it.

I have a table of points and I'm trying to generate a polygon around each point that estimates a maximum of 1000m driving distance from the point.

I'm using the current Ordnance Survey Vector Map District road layer loaded into postgis as multilinestrings, with Postgis 2.2 and Postgres 9.3.
All layers are using the OSGB 27700 projection.
The point data has a column which contains the nearest road node calculated from the source road data.

When I run the SQL for a single point to generate the list of individual roads using pgrouting driving distance it works perfectly:

SELECT seq,
       id1 AS node,
       id2 AS edge,
       cost,
       geom
  FROM pgr_drivingdistance(
       'SELECT gid AS id,
               source,
               target,
               st_length(geom) as cost
          FROM public.os_vd_2014_road',
       1529980, -- The node
       1000,    -- The driving distance
       false, false
                          ) as di
  JOIN os_vd_2014_road_nodes AS pt
  ON di.id1 = pt.id;

(Since this was essentially Anita Graser's code from her examples I expected nothing else).

However, I want to create a general case, ideally returning a polygon displaying the limits of each route but at the least a single geometry for each point.

To do this I created a function that would hopefully do that for a given point, and then call it in a simple SELECT statement:

CREATE OR REPLACE FUNCTION vdsl_hull (integer, integer)
RETURNS geometry AS $$

SELECT St_Collect(geom) AS the_geom
  FROM pgr_drivingdistance(
       'SELECT gid AS id,
               source,
               target,
               st_length(geom) as cost
          FROM public.os_vd_2014_road',
       $1,    -- The node
       $2,    -- The driving distance
       false, false
                          ) as di
  JOIN os_vd_2014_road_nodes AS pt
  ON di.id1 = pt.id;
$$ LANGUAGE SQL;

SELECT bt_pcp_list.pcp_id,
       bt_pcp_list.id AS gid,
       vdsl_hull(bt_pcp_list.nearest_node, 1000) AS the_geom
FROM   public.bt_pcp_list
WHERE  district_code = 'SMBRD';

This "appears" to work in that it generates output for each point, a valid "gid" value, and a column labelled "the_geom" but the geometry column is empty. The next stage was to surround the ST_Collect with ST_ConcaveHull() but clearly that's a stage too far for now.

I'm assuming the error in my thinking is something to do with the resulting geometry from the ST_Collect but could be completely wrong.

Any help would be very welcome!


texts: Finding if point exists in shapefile using shapely performance?

texts:
I have a data frame that consists of 2 million rows (lat, lon) and I am trying to figure out if each of these points are in a boundary defined by a shapefile or not (intact I am looking at multiple shapefiles). Here is the code, I have used to run this, but it is taking me around 30+ hours to run it. I am positive that there might be a better/faster way to do this and I am doing something incorrect.
I have also tried the approaches suggested here but to no improvement in performance (R-Tree) R-Tree Spatial Index, Aggregating Points

Here is the code that I am currently using:
To upload the shapefiles, I'm using geopandas, shapley, and pandas. I have downloaded the shapefiles from:
www.zillow.com/howto/api/neighborhood-boundaries.htm,
geo.nyu.edu/catalog/nyu_2451_34154, and
geo.nyu.edu/catalog/nyu_2451_34496

import fiona as fn
import shapely as sh
import geopandas as gpd
from shapely.geometry import Point
import pandas as pd

nyc_boroughs_df = gpd.read_file('./data/shapefiles/NYshapefiles/nyc_boroughs/nybb-geo_json.json')
nyc_neighborhood_df = gpd.read_file('./data/shapefiles/NYshapefiles/ZillowNeighborhoods-NY/ZillowNeighborhoods-NY.shp')
nj_neighborhood_df = gpd.read_file('./data/shapefiles/NJshapefiles/ZillowNeighborhoods-NJ/ZillowNeighborhoods-NJ.shp')
ct_neighborhood_df = gpd.read_file('./data/shapefiles/CTshapefiles/ZillowNeighborhoods-CT/ZillowNeighborhoods-CT.shp')
nyc_counties_df = gpd.read_file('./data/shapefiles/NYshapefiles/nyc_metro_counties/nyu_2451_34496-geojson.json')

def get_boro_city_county_name(row):
lat = row.lat
lon = row.lon
bn_series = nyc_boroughs_df[nyc_boroughs_df.geometry.contains(Point(lon, lat))].BoroName
if len(bn_series) == 0:
    bn_series = nyc_neighborhood_df[nyc_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
    if len(bn_series) == 0:
        bn_series = nj_neighborhood_df[nj_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
        if len(bn_series) == 0:
            bn_series = ct_neighborhood_df[ct_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
            if len(bn_series) == 0:
                bn_series = nyc_counties_df[nyc_counties_df.geometry.contains(Point(lon, lat))].name

if len(bn_series) > 0:
    return str(bn_series.iat[0])
else:
    return ''

df is a dataframe of 2 mil rows (columns: lat, lon) and when I apply the function on each of the rows, it takes around 30 hours to run.

df['name'] = df.apply(lambda x: get_boro_city_county_name, axis=1)


texts: Loading shapefiles from a read-only drive

texts:
Never had a problem with that, but after updating to ArcGIS 10.3 I cannot load shapefiles from a read-only drive. The error message says:

Error opening feature class
FDO error: -1 [nnn.shp]
FDO error: -1 [nnn.shp]

After clicking OK the following warning occurs:

Could not add the specified data object to the map.
FDO erro: -1 [nnn.shp]

Still I am able to open the shapefiles using other GIS tools. Probably ArcGIS 10.3 is trying to write a lock file into the folder. Everything worked fine in ArcGIS 10.0.


texts: What do ArcGIS regression analysis warnings 000853, 000916 & 000981 mean?

texts:
I'm trying to do regression analysis using a shapefile which contains dependent variable (NDVI) and explanatory variable (prec) in ArcGIS.
First, I run GWR analysis and then run spatial autocorrelation Moran I on residual of GWR result.
Finally, I wanted to see summary statistics on result window in ArcGIS but it just gives me 3 warnings related with analysis.
Do you have any ideas about why these warnings appear?
The shapefile has already been projected into Geographical lat/long which makes me wonder whether the result of regression is correct or not.


texts:
I have trouble with changing the image in layer control
I tried this: https://jsfiddle.net/v95h1kx4/

My browser has no problem with this but some people have top of the control blue and bottom of the control white - without the actual image. The weird thing is that they even use the same browser (Firefox, Chrome). It shows no errors for them too.


texts: Is a SQL String Aggregate possible in arcgis?

texts:
I have created a feature class from the results of a 1:many join between a table and a polygon layer in ArcGIS 10.3. The polygon layer contains land parcels and the table contains names. Where a 1:many join was present, ArcGIS created another row so for example if "John" and "Robert" both were owners, two rows would be created- one for John one for Robert.

I would now like to aggregate the results into a single field, so that I again only have 1 geometry containing a field with comma separated values for names e.g.

Parcelid: "10"
names: "John,Robert"


texts: Complete matrix from Generate Near Table

texts:
I am trying to make a matrix with all the distances from 4000 objects to 1000 other objects. I used the Generate Near Table and deselected "find closest option only". The output is table with 4 million rows. However I need it in excel in a 4000x1000 matrix.

Does anyone have an option to do this quick? Without creating 4 different excel tabs (maximum number of rows is 1 million) and using something like a INDEX MATCH function (this takes 24 hours to calculate per tab).

If you only need the closest feature, you can simply join the table. However, that doesn't work in this case.

Can anyone help me?

I have ArcGIS 10.1 and 10.2 and the object sets are both point polygons. I have spatial analyst.


texts: Reducing length of division and subdivision tics on ArcGIS Pro scale bar?

texts:
Using ArcGIS Pro 1.3.1 I have created a scale bar using the Scale Line 1 Metric template with a few modifications like:

Show one division before zero
Divisons: 11
Subdivisions: 10
Label Position: Below center
Label: KILOMETRES

It appears as below:

How can I make the division and subdivision tics (circled above) be only half as long as they currently are?


texts: ArcGIS Android - ArcGISFeatureLayer error: EsriServiceException: Internal Server Error

texts:
I am developping an Android application using tse ArcGIS SDK for Android. I am facing a problem when querying a featurelayer: the featurelayer I am using is hosted in our ArcGIS Online, and this is code I am using for the query:

public final String mNonValidatedURL="http://LINK-TO-SERVICE/FeatureServer/0";
QueryParameters qParam = new QueryParameters();
            qParam.setReturnGeometry(true);
            qParam.setOutSpatialReference(wm);
            qParam.setWhere("Validation='Non-validé'");
            qParam.setReturnGeometry(true);
            qParam.setOutFields(new String[]{"*"});
            Map<String, Order> obf = new HashMap<String, Order>();
            obf.put("Date", Order.ASC);
            qParam.setOrderByFields(obf);
            QueryTask qTask = new QueryTask(mNonValidatedURL);
try
            {
                FeatureResult results = qTask.execute(qParam);
                return results;
            }
            catch (Exception e)
            {
                Log.e("freaking exception", e.toString());
                e.printStackTrace();
                return null;
            }

sometimes the query is executed without problems but most of the time I am catching this error:


texts: Concatenating date fields with Field Calculator ArcPy

texts:
I am attempting to calculate a field based on date and time of two other fields.  How should my expression read in ArcPy?  Both fields are date type with the values of the ORIGDTDATE reading 11/18/2016 while ORIGDTTIME reading 9:52:00 AM  I would like for my new output for ORIGDTDATE to equal 11/18/2016 9:52:00 AM

import arcpy
calcExpression =  [!ORIGDTDATE!] + [!ORIGDTTIME!]
arcpy.CalculateField_management(fc, "ORIGDTDATE",ORIGDTDATE + !ORIGDTTIME!,"PYTHON_9.3")


texts: Matching Decedal Average NDVI and Worldclim Data

texts: Ive never used NDVI (or EVI) datasets in R. Does anyone have any resource I could use to access a decadal average NDVI dataset for the globe. It would be preferable if this dataset was also comparable in dimensions and resolution to the Worldclim bioclim variables (2.5 arc mins) so that the rasters can be easily stacked. Any directions would help here.


texts: Feature Class Update Trigger with Archiving

texts:
We are deploying Collector in our institution and for offline work on unversioned data we need to enable archiving.

The feature class needs some updates when the user enters new records using the collector or webmap.

So to make this update we put a trigger on the archiving-enabled feature class to update the view of that feature class, as the documentation says.

The problem is that when this trigger is fired, it appears that the record was not fully created in the feature class, and the result is that the update made by the trigger is not applied to the register with the gdb_to_date attribute value 12/31/9999 but in the previous record.

In an answer at GeoNet to SQL Trigger on table causes CODE 500 Error, someone said that when an insert is made into an archived feature class, an insert is made and then an update, and the trigger is firing on the insert.

Our database is PostgreSQL.

Below the records that represent the problem described, note that the attribute "situacao" is updated with the value "ATIVO" by the trigger and then the final record is created (gdb_to_date 12/31/9999) without the value "situacao" = "ATIVO"

In my view the update trigger should be fired after this last record gdb_to_date 12/31/9999, but that's not what happens.

Nor do we understand why the update is deleted from the final registry.


texts: Are geotiffs with the embedded georeferencing information supported?

texts:
I am able to load raster data into my postrgreSQL data base however when I display in Quantum GIS the pixel height and width of the image is being translated into real world lat/lon coordinates, as a result my image appears 1000 degrees wide and 600 degrees high which is obviously wrong. The Image I am loading is a geotiff ie the coordinate header is embedded in the file rather than a world file.

I am using the following combination of database and software:

Quantum GIS 1.8.0
Postgis 2.0
PostgreSQL 9.1

Can anyone elaborate or point me in the direction of an explanation of how to get the image to load correctly?

Update:

This is now working for me as I have tried loading using a tiff file which includes the separate world file (.tfw). So I take it geotiffs with the embedded georeferencing information are not supported?


texts: Changing marker shape based on a rule or column value in QGIS?

texts:
I am using QGIS 2.14.14. I am styling point markers in a layer that I imported using "Add Delimited Text Layer". I want to style the markers based on three properties:

Color of marker relates to study name
Size of marker relates to number of soil samples at that location

(I have successfully completed the first two requirements, especially with the help of this post:
How to style points in QGIS based on multiple attributes).

I used a categorized style and made a color gradient to create the different colors for the studies, and I used the "size assistant" to change the size of the marker.

3. Shape of marker denotes if the soil samples were taken in "topsoil" or "other".

Specifically, I would like the shape to be a square if it's taken from "topsoil" and a circle if it's taken from "other".

I can not find a way to change the shape based on a condition or a certain value in my attributes table. I have tried making a column that's a whole number field and binary (1 for topsoil, 2 for other). I have also tried making a text field ("Topsoil" for topsoil and "Other" for other).

Does anyone know of a way to successfully combine these three formatting conditions? Specifically, does anyone know how to conditionally format a marker shape based on some attribute?


texts: Append Annotation within a ArcGIS geodatabase

texts: I have a large data set of sanitary sewer pipes. I am making a map book with data driven pages and need to move the labels around so that they are not overlapping the edges of the page and the manhole and catch basin labels are legible as well. I have converted the labels to annotation in the database and can easily edit the labels, moving them. I noticed that not all the features were labeled, so I check the unplaced annotations. Everything is supposedly placed. I realized that I had the connect features option check in Placement Properties > Label Density prior to converting to annotation. I believe this is why many of the labels are not showing up. Is there a way to update or fix this? I have about half the labels moved already and don't want to start over.


texts: Assigning Values to Consecutive Line Features

texts:
I'm trying to assign limiting capacities to pipes that are in series and also branch off of each other. If the current pipe's capacity is > the previous pipes capacity, then the current pipe's limiting capacity = the previous pipes capacity.

In the screenshot above, the pipe circled in red would be the limiting capacity for the pipe that is to the right of it, however it wouldn't be the case for the pipe below it.

As of right now I've been able to select these pipe branches based on their connection points (Manhole IDs) using the following code:

#Selects pipes that are intersecting trunks
arcpy.Select_analysis (links, "Trunks", """ "TRUNK" = 'YES'""")
arcpy.SelectLayerByLocation_management(links, "INTERSECT", "Trunks")
arcpy.SelectLayerByAttribute_management(links, "REMOVE_FROM_SELECTION", """ "TRUNK" = 'YES'""")

#Creating a dictionary that stores the US (Inlet) and DS (Outlet) manholes and the struct_ID
system = dict()

with arcpy.da.SearchCursor(path, ["INLETNODE", "OUTLETNODE", "NAME"]) as cursor:
for row in cursor:
    if row[1] not in system:
        system[row[1]] = [],[]
    system[row[1]][0].append(row[0])
    system[row[1]][1].append(row[2])

#Finds the pipes upstream of the selected pipe and adds them into the us_pipe_list list
us_pipe_list = []

def find_us(inletnode):
    if inletnode in system:
        us_pipes = system[inletnode][1] #Struct_ID of pipes
        for us_pipe in us_pipes:
            if us_pipe in us_pipe_list: #If the Struct_ID of a pipe is in the list, move on
                return
            else:                       #If the Struct_ID isn't in the list:
                us_pipe_list.append(str(us_pipe)) #add it in,
                us_mh = system[inletnode][0] #and US manhole of pipe
        for mh in us_mh:
            find_us(mh)

pipe_layer = arcpy.mapping.Layer(links)
with arcpy.da.SearchCursor(pipe_layer, "INLETNODE") as cursor:
    for row in cursor:
        find_us(row[0])

where_clause = "NAME in ("+str(us_pipe_list)[1:-1]+")"
arcpy.SelectLayerByAttribute_management(links, "ADD_TO_SELECTION", where_clause)

What I've thought of so far is to add in the Pipe Capacity field into the system dictionary:

with arcpy.da.SearchCursor(path, ["INLETNODE", "OUTLETNODE", "NAME", "PIPE_CAPACITY"]) as cursor:
for row in cursor:
    if row[1] not in system:
        system[row[1]] = [],[],[]
    system[row[1]][0].append(row[0])
    system[row[1]][1].append(row[2])
    system[row[1]][2].append(row[3])

And then call on these Capacities through the same find_us function that's used to add the pipes that are to be selected to the list (us_pipe_list) that is called on when using the actual Select by Attribute tool at the end. Below is the snippet of code that I'd add in to the find_us function after the us_pipe for loop:

lim_caps = system[inletnode][2] #Capacity of pipes
    for cap in lim_caps:
        if cap > cap_1:
            cap = cap_1
        else:
            cap = cap
        cap_1 = cap

With all that being said, I'm not sure how to actually populate the limiting capacity field after comparing the previous pipe's capacity (cap_1) and the current pipe's capacity (cap) to see if cap > cap_1.


texts: Select by Attributes not selecting ALL larger numbers greater than 75000 in ArcMap?

texts:
I was trying to select all values over 75000 in Select by Attributes.

Surely I thought doing so would do just that, but it's essentially selecting only values that begin with 75 and 8 and 9 (even includes a low value like 900) and leaving out values like 100000, 200000, and 700000. I entered this:

SDE_DATE.DBO.Parcels.Total_AV > '75000'

Do you know how I can select all values greater than 75,000?


texts:
I would like to import an esri Personal Geodatabase (.mdb) to PostGIS.
I have tried to follow the blog post:

How to import esri Personal Geodatabase (.mdb) to PostGIS on Linux

I installed the GDAL tool suite (using gdal-110-1600-x64-core.msi), but whenever I try to run ogr2ogr, I get the following error:

ogr2ogr -f "PostgreSQL" PG:"host=localhost user=penusr dbname=my_pg_db password=my_pg_db_pass port=5432" C:\data\PROVINCES_2013.mdb -a_srs EPSG:26986
ERROR 1: Unable to initialize ODBC connection to DSN for DRIVER=Microsoft Access Driver (*.mdb);DBQ=C:\data\PROVINCES_20
13.mdb,
[Microsoft][ODBC Driver Manager] Data source name not found and no default driver specified
FAILURE:
Unable to open datasource C:\data\PROVINCES_2013.mdb' with the following drivers.
  -> ESRI Shapefile
  -> PostgreSQL
... list of many more.

Does this error typically indicate a problem with connecting to my Postgre DB, or something else? Am I missing some settings? I am fairly sure that the dbname and password are correct.

If I go to the ODBC Data Source Administrator window, I can see that under the Drivers tab PostgreSQL ANSI(x64) and PostgreSQL Unicode(x64) are listed.