texts: Layer visibility using maxscale minscale not working

texts:
I'm using OpenLayers and GeoServer. I have a series of layers that have a scale parameter like: 1:50000-1:25000 where 50000 is the max scale while 25000 is the min scale.

I have a map like in this example:

map = new OpenLayers.Map(
{
        allOverlays: false,
        projection: "EPSG:32632",
        maxExtent: bounds,
        scales:[500000,350000,250000,100000,25000,20000,15000,10000,5000],
        units: 'm'
});

And I have this layer WMS:

var layer = new OpenLayers.Layer.WMS("layer1",wms_address,
        {
            layers:"ws:layer1",
            transparent:true
        },
        {
        maxScale:50000,
        minScale:25000
    });

and I would expect that the layer within that scale it doesn't show up in the map. I need this kind of behaviour. However if I set the scale 1:25000, the layer is present inside the map (I saw it with Firebug) but it is not showed up in the map!

I also tried with max and min Resolution but nothing, the layer is in the map but I cannot see it.

What I'm doing wrong?

Please help me, I looked everywhere in every forum I don't know how to solve this problem!


texts: I'm trying to map the volume of people travelling from various cities to popular mountain towns in British Columbia. I have a CSV containing the starting point, ending point and number of travellers. There must be a way to show the volume of travel (thickness of line, or colour coded) from place to place using city layer polygons and a road network layer.


texts: exporting a featureclass to shapefile with selected fields only

texts:
Possible Duplicate:
  How to control which Fields are exported in ArcObjects? (via Field Mapping)

I want to export a featureclass to a shapefile.The only condition i want to impose is that shapefile should contain only selected field from featureclass.Somebody plz help out.


texts: How to convert a polygon to a polyline using the sf package in R?

texts: I have some simple polygons that I would like to convert to polylines using the sf package in R. I have found st_polygonize to convert polylines to polygons but I cannot but the reserve operation. Some documentation suggested that st_segmentize would work but I haven't found it useful.


texts: Getting XY map coordinates on mouse clicks as input variables into a python script

texts:
This is my first time coding in python, so bear with me. I wrote a script to calculate an analytical solution to find the orientation of a plane based on three points. I'm using ArcGIS 10.1.

# Input points
x1 = float(arcpy.GetParameterAsText(1)) # Point 1
y1 = float(arcpy.GetParameterAsText(2))
z1 = float(arcpy.GetParameterAsText(3))
x2 = float(arcpy.GetParameterAsText(4)) # Point 2
y2 = float(arcpy.GetParameterAsText(5))
z2 = float(arcpy.GetParameterAsText(6))
x3 = float(arcpy.GetParameterAsText(7)) # Point 3
y3 = float(arcpy.GetParameterAsText(8))
z3 = float(arcpy.GetParameterAsText(9))

# Plane equation
A = (y1*z2) + (z1*y3) + (y2*z3) - (z2*y3) - (z3*y1) - (z1*y2)
B = (z2*x3) + (z3*x1) + (z1*x2) - (x1*z2) - (z1*x3) - (x2*z3)
C = (x1*y2) + (y1*x3) + (x2*y3) - (y2*x3) - (y3*x1) - (y1*x2)
D = (z1*y2*x3) + (z2*y3*x1) + (z3*y1*x2) - (x1*y2*z3) - (y1*z2*x3) - (z1*x2*y3)
E = np.sqrt(A**2 + B**2 + C**2)

AZ = np.arctan(A/B) * (180/np.pi) # Preliminary azimuth for dip direction
dip = np.arcsin(-np.cos((np.pi/2) + np.arccos(C/E))) * (180/np.pi) # True dip

# Test for quadrant
alpha = A/E
beta = B/E

# Place dip-direction in proper quadrant
if alpha > 0 and beta > 0:
  DD = AZ
elif alpha > 0 and beta < 0:
  DD = 180 + AZ
elif alpha < 0 and beta < 0:
  DD = 180 + AZ
elif alpha < 0 and beta > 0:
  DD = 360 + AZ

# Right-hand rule strike
if DD - 90 < 0:
  RHR = 360 + (DD - 90)
else:
  RHR = DD - 90

# Centroid coordinates
Cx = (x1+x2+x3)/3
Cy = (y1+y2+y3)/3
Cz = (z1+z2+z3)/3

# Insert centroid point
cursor = arcpy.InsertCursor(arcpy.GetParameterAsText(0)) # Feature class must already exist
row = cursor.newRow()
row.Shape = arcpy.Point(Cx,Cy,Cz)

# Populate table
row.X = Cx
row.Y = Cy
row.Z = Cz
row.DD = DD
row.dip = dip
row.RHR = RHR
cursor.insertRow(row)
del cursor, row

# Refresh MXD
arcpy.RefreshTOC()
arcpy.RefreshActiveView()

This works, except that the user has to hover each point and manually enter the coordinates. Ideally, I would like:
(1) the user initially define the target feature layer in a pop-up menu, and hit OK;
(2) click on the map to extract the first XY map coordinates;
(3) pop-up queries the user for the first Z, hit OK;
(4) click on the map to extract the second XY map coordinates;
(5) pop-up queries the user for the second Z, hit OK;
(6) click on the map to extract the last XY map coordinates;
(7) pop-up queries the user for the last Z, hit OK; then
(8) the script runs to completion.

I found the ESRI tool tutorial that allows the map XY to be displayed on a message box after a mouse-click:

import arcpy
import pythonaddins

class MC(object):
    """Implementation for MouseClickSample_addin.MCtool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
    def onMouseDownMap(self, x, y, button, shift):
        message = "Your mouse clicked:" + str(x) + ", " + str(y)
        pythonaddins.MessageBox(message, "My Coordinates")

I cannot figure out how this Python Add-in can pass the XY variables as input to my script.


texts: Most optimal algorithm for contour correctness?

texts:
Am creating an app which correctly calculates contours from a DEM file or any height-map given to it. The correctness of the contour holding the most importance, since a person viewing the contours shouldn't end up taking wrong planning decisions and maybe losing their life by the route they choose after looking at the map.

I have already tried using GDAL 1.7 source code, compiled it and added some features of my own to it for some other objective. Now I want to add-to/modify the source to be able to use the contour generation function of GDAL in C++.

But before proceeding, I wanted to know:
1. If GDAL offers the best free contour generation algorithm available or is there any other library which I can link to my code and feed my DEM map to it and receive contour vector lines (or go thru research papers and implement the contour generation algorithm myself in C++)? If GDAL, then why is it the best?
2. A generated contour will basically be points connected with straight lines. Is it more correct for a contour to be depicted with connected straight lines or will it give a more correct result to use a spline to smoothen out the lines?

I won't be interpolating any contours or the DEM because I've read on this website that interpolation leads to incorrect results.


texts: Converting multiple GPX to SHP files using ModelBuilder in ArcMap

texts:
I usually use the "GPX To Features" tool to import tracks from my Garmin GPS device into ArcMap. Now I have a folder which contains around 50 GPX files which I want to convert to SHP files in one go. Using the ModelBuilder I created this model:

The problem is, that I cannot connect the green "File" oval (which should contain the current file name) as an input with the "GPX To Features" tool. The only available type of connection is "Precondition". What am I doing wrong?


texts: LISA (local indicators of spatial association) p values in R versus in Geoda

texts:
When calculating LISA statistics in GeoDa versus using the localmoran() function in spdep, the I values are nearly the same, but the p values can be quite different.

I know the two programs calculate significance in different ways, but I'm not sure I can clearly explain how. Can someone suggest a good way to explain to a mixed audience (to those with advanced and those with basic statistical knowledge, and to someone like me who's kind of in the middle) why the two results are different?

As the significance results are so different, are there reasons to prefer one method of calculation over the other?

Is there a way to make the results from the two programs more similar?


texts:
One might expect that the "Color" CAD field would control this, but the driver doesn't allow it to be set.

schema = {'Layer': 'str', 'Color': 'int'}

with fiona.open(path, 'w', 'DXF', schema=schema) as f:
    f.writerecords(records)

Gives:

ERROR: failed to export vectors: DXF layer does not support arbitrary field creation, field 'Color' not created.

The driver documentation suggests that it attempts to interpret OGR feature styles, but I see no mention of these in the Fiona docs.  Is it possible to set arbitrary colors with Fiona?


texts: Rendering elevation data of planet_osm_line in GeoServer?

texts:
I have created a postgres data from osm data of sachsen and from NASA srtm data extracted elevation data of the same region. I have multilinestring with varying elevation data for a single road, so is it possible to display elevation by using geoserver.

From what I know there is an option of additional dimension for elevation that can be enable in geoserver with a particular column set as the attribute that corresponds to elevation in the database.
 But since my lines are multilines each road can have varying elevation points can I display varying elevation on the line.

Can this be set for multiline?


texts: Map doesn't appear with custom CartoDB query

texts:
I'm trying to use CartoSQL from CartoDB.js to create a layer. Here is the SQL query I am trying to make:

SELECT
  rpa_nj_hsip_hospitals_compressed.flood,
  rpa_nj_hsip_hospitals_compressed.the_geom
FROM
  rpa_nj_hsip_hospitals_compressed

All the columns are in the table –– I checked. This is the query string and you should see the same error I am seeing.

What am I doing wrong?


texts: Reducing size of KML file (after converting from Shapefile)

texts:
I have a shapefile of size 224 MB. I am using ArcGIS to convert it to a KML file. The process involves converting the shapefile to a feature layer and then converting the layer to KML.

An alternative method that I have tried is using ogr2ogr.

In both the cases, the KML file produced is around 500 MBs. I need it to be less than 250 MBs, so I can import it as a Google Fusion table.

I was informed that one of the best ways to cut size is to reduce precision of the coordinates to 5 decimal places (from its current 10 or so).

How do I go about changing the decimal places using either ogr2ogr or ArcGIS? Is there some other software that may help in this and also, are there any other strategies I may use to reduce the size of the KML file?

I tried using Notepad++ to open the KML file and then use regular expressions to change coordinate precision but for some reason this method is not working. If I make the change and try and save the file, all the data gets deleted and the KML file size changes to 0 bytes.

An alternative strategy that I was thinking of was to split my shapefiles into two or three files and convert them to KML separately, upload them to Google fusion tables and then append those. However, I am not sure of two things.

How to split shapefiles into two or three smaller shapefiles?
Whether I can easily append fusion tables together to create a larger fusion table?


texts: How to wire ESRI.ArcGIS.Desktop.AddIns.Extension to a custom AddIn Button click event?

texts:
I am developing custom ArcMap toolbar with various AddIn tools, button and extension. I see the sample how to access the ActiveView events in the extension Initialize() method.

 private void Initialize()
    {
      // If the extension hasn't been started yet or it's been turned off, bail
      if (s_extension == null || this.State != ExtensionState.Enabled)
        return;

      // Reset event handlers
      IActiveViewEvents_Event avEvent = ArcMap.Document.FocusMap as IActiveViewEvents_Event;
      avEvent.ItemAdded += AvEvent_ItemAdded;
      avEvent.ItemDeleted += AvEvent_ItemAdded;
      avEvent.SelectionChanged += UpdateSelCountDockWin;
      avEvent.ContentsChanged += avEvent_ContentsChanged;
             ....
              // ? Possible to access AddIn button even here ?
              myCustomAddInToolbar.myCustomAddInButton.RaiseCustomEvent += HandleCustomEvent;

How can I wire my custom AddIn Button? I don't know how to access the AddIn button events from here.
Is there any ArcMap(ArcObjects) Interface to access AddIn Toolbar/Button Events?
Thanks


texts: Proj4js - no transformation for epsg2056

texts:
I have a script which converts local coordinates to WGS84.
The script works fine if source data and source CS is epsg:21781. The CS for the data has changed now to epsg:2056. When I use epsg:2056 as source the point coordinate is not transformed at all and the transformation output equals transformation input. Here is the script where the transformation takes place:

var src = new Proj4js.Proj("EPSG:21781"); //<-- EPSG:2056 does not transform
var dst = new Proj4js.Proj("WGS84");
var point = new Proj4js.Point(center.lon, center.lat);
Proj4js.transform(src, dst, point)

I don't know anything about Proj4js and I'm just stuck here.


texts: Display google map with boundary lines for each address

texts:
I am working on Google MAP API integration and I have tried to display the map with boundary lines for each address, such as Zillow.com.

Eg: http://www.zillow.com/homes/for_sale/priced_sort/25.945776,-80.147139,25.9432,-80.152374_rect/17_zm/

Zillow gives a neighborhood api/db file with boundaries and I'm not able to make sense of it exactly. How to get map view like their view with the DB data?

I have addresses with latitude/longitude. Using these addresses and latitude/longitude can I draw the boundary line for each address? If the address specifies a building then it needs to display it with its boundary line. Is there any API from google MAP or third party?

Please give me suggestion on this - I went through google search results and related articles on this and couldn't get any result.


texts:
I am trying to process a large number of HRRR files, about six months' worth of them, which is close to 4059 files. I am working in R, and this is the function I have, that I am trying to vectorize over all the files in list.files('path to HRRR files, full.names = T):

    hrrr_processing_fn <- function(hrrr_file,vars){
  file_date <- str_extract(tail(unlist(str_split(hrrr_file,"/")),1),"[^\\D$]{8}\\.[^\\D$]{2}")

  message(hrrr_file)

  # 1) Identify the layer number of each of the variables
  index.file <- varfiles[grepl(file_date,varfiles)]

  hrrr.idx <- fread(index.file, sep = ":", header = F)
  hrrr.idx[,V5 := str_replace_all(V5, " ","_")]
  var_num <- hrrr.idx[V5 %in% c("surface", "2_m_above_ground") & V4 %in% vars,V1]
  colnames <- hrrr.idx[var_num,paste0(V4,"_",V5)]

  # 2) Read the grib file as a RasterStack
  hrrr_stack <- stack(hrrr_file, bands = var_num)

  # 3) Reproject the rasterstack to the required CRS
  hrrr_stack_crs <- projectRaster(hrrr_stack, crs = crs_latlon$proj4string, method = "ngb")

  # 5) Mask the rasterstack to the shapefile
  hrrr_stack_masked <- mask(hrrr_stack_crs, mask = suffolk_county)

  # 6) Extracts the centroid and data of the masked raster as a data.table
  hrrr_mask_dt <- as.data.frame(hrrr_stack_masked,na.rm = T,xy = T)
  setDT(hrrr_mask_dt)
  setnames(hrrr_mask_dt,c("lon","lat",colnames))
  file_date <- str_replace(file_date, "\\.","_")
  hrrr_mask_dt[,date := ymd_h(file_date, tz = "America/New_York")]

}

where crs_latlon = Coordinate Reference System:
  EPSG: 4326
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

and vars = 'TMP'

While most files are being processed, there are some files which throw an error message, such as this:
Error in .local(.Object, ...) :
  /home/my_username/path_to_files/20160806.22.grib2 is a grib file, but no raster dataset was successfully identified.

This is not the only error message I get, though. Some files are just not readable. I have ensured that the files are present and downloaded correctly, because I got a '200' message when I downloaded the file using download.file().

These errors halt the function from running. Given the large number of files, I do not want to stop at each bad file and start over again. I want to write a way into my function to not halt the processing, but log the files that are 'bad' in a different file, and keep processing the other files. One way I could think of was to use a conditional statement examining if the file was readable and read without problems, in which case to go ahead with the processing, or if the file could not be read, which would lead to logging the file name, but I am unsure how to operationalise it. I have tried the assertive package function is_executable_file, which does recognise the grib2 file in question to be an executable file, but the file is nevertheless not read by raster::stack. I have a vague idea that I need to use tryCatch for this, but my attempts so far have failed. Any help would be greatly appreciated.  Thanks!


texts: ArcMap Maps Show in Layout View but Not Data View

texts:
I can see all my map layers in Layout view but get a blank white screen in Data View. Every help resource I searched seems to show thousands of questions and answers for the opposite problem.

I have tried zoom to layer and zoom to visible. I have tried setting display to wire frame and to always display. There are no red x's on any of the map layers and the source path's have not changed. These were all working fine yesterday.

This morning it was good, then I zoomed in and got a blank layer, zooming back out didn't return to the previous map. Pulling my hair out. Using military map deck made for CPOF, if that information is helpful to anyone.

Version 10.1
Problem exists on all saved and new mxds.
Problem exists at all scales.


texts: Rendering a QGIS WMS layer in Firefox hogs the machine

texts:
I am less than a novice beginner on GIS Systems and am groping for answers.

I followed instructions for WIN 7 Install at both
http://underdark.wordpress.com/2012/04/06/qgis-server-on-windows7-step-by-step/ and http://hub.qgis.org/projects/quantum-gis/wiki/QGIS_Server_Tutorial#Windows

All appears to work well.  I can retrieve a capabilities-XML.

Next, I add an WMS Layer to a shapefile in QGIS.  Again, no issues.

Then I attempt to display via Firefox 18.0.2 with address call

  [http://localhost]/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&SRS=EPSG:4326&REQUEST=GetMap&map=C:/OSGeo4W/apps/qgis/bin/VA_VOTES_BUSINESS.qgs&BBOX=-36.215710,6.817432,-9.700450,58.551008&WIDTH=500&HEIGHT=1000&LAYERS=VA_VOTES_BUSINESS&FORMAT=image/png

The entire operating system becomes unstable at the page call.  Firefox displays Connected to LocalHost; there is 100% CPU usage (Pentium Dual Core 2.2 GHz); and 80% memory usage (4 GB).  I have to power off to regain control.

I see nothing but errors in the Apache error log but of course I'm not familiar with the references.  Spent 10 hours+ trying to troubleshoot and this is my last resort.  Any direction would be most helpful and appreciated.

QGIS 1.8.0 Lisboa Desktop


texts: GDAL - Command Line vs. Spyder (Python) Functionality

texts:
I am trying to retrieve information from a GeoTiff using the GDAL API in Python (Spyder).  Within Anaconda I have my base environment which contains gdal v2.3.3, but I also created a virtual environment with gdal v3.0.2 since this version of GDAL can read GeoTiffs that can be classified as "BigTiffs" (> 4GB).
When I run 'gdalinfo filename.tif' from the Anaconda terminal within my virtual environment, I get coordinate system information output that starts with:
Coordinate System is: PROJCRS["WGS 84 / UTM zone 16N", BASEGEOGCRS["WGS 84", DATUM["World Geodetic
System 1984", etc..

I can also obtain this same data when I run 'print(gdal.Info("filename.tif"))' within Spyder which is also configured for my virtual environment:
from osgeo import gdal, osr, ogr
print(gdal.VersionInfo())  # outputs 3000200 (i.e. v3.0.2)
print(gdal.Info('filename.tif'))  # outputs same data as Anaconda terminal

However, when I attempt to obtain GeoTiff information using 'gdal.Open("filename.tif")' within Spyder and invoke other GDAL functions I obtain GeoTiff data that does not appear compatible with gdal v3.0.2:
dem = gdal.Open('filename.tif')
print(dem.GetProjection())  # outputs PROJCS["WGS 84 / UTM zone 16N",GEOGCS["WGS 84", etc...

For the "GetProjection()" function, for example, I need the output to show "PROJCRS" and not "PROJCS".  I need the GDAL commands within Spyder to be compatible with v3.0.2, but they seem to be returning outputs compatible with v2.3.3.  Is there something else I may be missing when invoking the GDAL API?


texts: Using Calculate Field to fix hyperlinks in ArcPy gives Cannot use VB for services?

texts:
Getting an error while trying to script it, specifically "Cannot use VB for services."  Snippet and VB code block as follows:

Raw block (works within ArcMap): "\ftp\raw\YYYYMMDD\Files\Charles" +[NAME]+".ntf"

Snippet:

    date = '"20160207"'

    arcpy.CalculateField_management(outShp, "Link", '"\\\\ftp\\raw\\"' + date + '"files"' + "[NAME]"+ '".shp"', "VB", "")

Not sure what the issue is, somehow need to be able to insert the date and have it syntactically be correct (in either VB or python).

I didn't explain very well I have a column that contains a portion that is part of the hyperlink, so need to concatenate A column to the end of B column (with the partial path) to complete the link


texts: How to create a rectangle polygon with set lengths in QGIS?

texts:
I have a search area i want to outline. and i need it to be 1km by 2km. is it possible to achieve this in QGIS.

I can create a rectangle polygon but only guess the lengths.


texts:
I'm trying to run SELECT in order to extract specific features under "CLASS_NAME", but continue getting "ERROR 000210: Cannot create output".

I've searched though other questions with the same error but to no avail.  https://gis.stackexchange.com/search?q=ERROR+000210

from arcpy import env
env.workspace = "C:\input"
arcpy.Select_analysis("oldpoly.shp", "C:\output\newpoly.shp", '"CLASS_NAME" = item1,item2,item3')

Still getting to know Python and learning how to troubleshoot.


texts: Import csv value to an empty field in a point feature class

texts:
I have a .csv with three colums and a point feature class with an empty field that i would like to fill with a value of the csv. In the .csv i have 15 rows with a code. In the feature class just one row and one of the fields has the same code as one of the codes in the csv.

I want to read the feature class field with code, then read the csv and when it find  the same code, change the colum in the csv and read the value corresponding to it. And that value which correspond to the coincident code value insert it in the empty field in the point feature class. By now i was trying to use the SearchCursor and csv reader but it seems it doesn't work but also I get no error.

Does anyone know how to do that?

import csv
import arcpy
fc = r'C:\Users\Alumno\Desktop\IGP\IGP.gdb\valor'
with arcpy.da.SearchCursor(fc,['Com__Forest']) as cursor:
...   for row in cursor:
...      comarca = 'row [0]'
...      with open (r'C:\Users\Alumno\Desktop\IGP\FWI.csv','rb')as f:
...          reader = csv.reader  (f,delimiter = ';')
...          if row in reader == comarca in row:
...                 indiceCom = row[2]
...                 print row [2]

my problem is that it is the first time i am leading with python and i am terribly lost. From one side i arrange to compare the value on the csv with one difined value of "Comarca" import csv

comarca = '1'

with open (r'D:\ProyectoIGP_Inigo\FWI.csv','rb') as f:

... reader= csv.reader (f,delimiter = ';')

... for row in reader:

... if row[0] == comarca in row:

... indicecom = row[2]

... print indicecom ...
1

Now i need to change that defined value of "comarca = '1'" for and edit cursor or something that lets me introduce the value of row [2]

Any idea?


texts: Build Radio Button in QGIS Edit Widget Properties

texts:
My question is similar to Choosing multiple options from QGIS Check Box or Value Map for feature attribute?

I have 3 categories in a row, i want to make this selection in form of radio button so that selection can be made more than one.

Instead of using Value Map just like below, I want to display it using Radio Button.

Can you guys help me out?


texts:
I have a layer.

I want to transform the type of geometry from point to polygon while creating a buffer around points proportional to the population of cities.

I'm testing :

ALTER TABLE city
ALTER COLUMN points TYPE POLYGON USING ST_buffer ('points', pop * 100)

but that does not work...

Can anyone put me on the right track?


texts: Calculating field to equal adjacent polygon field when field = 0

texts:
I have a polygon Feature Class in ArcGIS Pro 2.5.0

In the image below, yellow is when the field value is 0, and blue is all other values.

I want to calculate all the 0 values to equal the adjacent polygon value of the same field.


texts:
I'm using NetTopologySuite to read a shapefile like below code:

        private FeatureDataTable GetFeatureTableFromShapefile(string shapeFilePath)
    {
        GeoAPI.GeometryServiceProvider.Instance = new NetTopologySuite.NtsGeometryServices();
        FeatureDataSet ds;
        using (var sf = new ShapeFile(shapeFilePath) {Encoding = Encoding.Default})
        {
            sf.Open();
            var ext = sf.GetExtents();
            ds = new FeatureDataSet();
            sf.ExecuteIntersectionQuery(ext, ds);
        }

        return ds.Tables[0];
    }

Now I want to merge all features (GeoMetry) and create only one feature from them.

Please note that I'm trying to create one feature and not something like multipolygon like below image:

For above image, i used this code, but it's not what i want, because there still a red line between to features that i tried to merge:

            var entities = new List<DbGeometry>(){ d1.Geom, d2.Geom};
        DbGeometry allBoundaries = null;
        for (var i = 0; i < entities.Count; i++)
        {
            if (i == 0) allBoundaries = entities[i].Boundary;
            else allBoundaries = allBoundaries.Union(entities[i].Boundary);
        }

So is there any way to accomplish this?


texts: ArcGIS Server Error: “Feature service requires a registered database”

texts: I'm attempting to publish a feature service to ArcGIS server. I'm running PostgreSQL and ArcSDE. When I try to publish the service I get the error "00090: Feature service requires a registered database". When I validate the database registration it appears to succeed, but I still can't publish the service. The features I am publishing are located inside of the geodatabase that I am registering. This geodatabase is located on the same server that ArcGIS Server resides. Am I missing something ?


texts: Plotting raster over openstreet map after reading in with rasterio

texts:
I have read in a simple raster with rasterio that I would like to plot over an OSM map (interactive if possible).
- My first hurdle is plotting raster with real world coordinates rather using default 0,0 origin.
- Second is finding library to plot interactively on (folium, cartopy, etc). I have failed to implement raster overlay on these

Here's the ipynb: https://github.com/arnawldo/geospatial-scratch/blob/master/rasterio.ipynb


texts: How to adjust KML files so that they display properly in Google Earth, ArcGIS Desktop & ArcGIS Online?

texts:
The students in one of our labs have been cranking out daily satellite derived Sea Surface Temperature and Chlorophyll maps for a few years. They are published for public consumption and we have been using them in our global visualization lab overlaid on Google Earth for years.

We have been making them available on a lab website in the browser overlaid on Google Maps and as a Google Earth plugin. Now we would like to embed them into an ArcGIS Online web app so they can be more easily embedded on other lab pages.

The problem is that ArcGIS Online appears to be shifting the PNG overlay images, but a colleague of mine who has ArcGIS Desktop says that the images are not shifted when he opens the KML files using them using Desktop.

We opened a ticket with ESRI support and the answer we got back was:

"The reason for this discrepancy comes from the difference in projections between Google Earth and ArcGIS Online. Google Earth uses a "Simple Cylindrical projection with a WGS84 datum for its imagery base" [see reference: https://support.google.com/earth/answer/148110?hl=en] while ArcGIS Online uses WGS 1984 Web Mercator (Auxiliary Sphere). While these projections are similar, they're not the same, and therefore you will not see the exact same results in both maps when displaying the same data." And that "Although you are correct in that the data does not align in AGOL identical to Google Earth, the data is still be displayed as-expected, so I won't be able to log a bug or enhancement request for this situation."

While I imagine we could try to alter our code so that the PNG file is based on a different projection, I'm thinking that will hose those users who are referencing them for their research as now they'll have a problem displaying the KMLs using Google Earth. Plus I'm confused why ArcGIS Desktop would display them properly and ArcGIS Online would shift them.

The PNG files that the KMLs reference are updated daily and are available here:
http://modata.ceoe.udel.edu/dev/mcimino/MODIS/GSFC/KMZ/
and
http://modata.ceoe.udel.edu/dev/dhaulsee/MODIS/GSFC/KMZ/

Any thoughts as to how we can resolve this so that the existing users of the data aren't disrupted but where ArcGIS Online web apps can be leveraged to facilitate expanded web visualizations of the data? It seems naive if ArcGIS Online claims to display KMLs, but only if they're re-projected specifically for use on ArcGIS Online only.


texts:
I have an script to transform DN into Reflectance of each band of several images.

The unique variable of the equation, apart of the Band itself, is the Divisor:

(B1 - 1) / 508
(B2 - 1) / 254
(B3 - 1) / 363
(B4 - 1) / 423

I wrote an script, copying the same code four times.

I know it can be improv it, using a function instead of repeating it, but I do not figure out how to do it.

Any idea?

import os, arcpy, glob
from arcpy import  env, sa
check = arcpy.CheckOutExtension("Spatial")
print check

arcpy.env.workspace = r'F:\Mosaic 2000'
ws = env.workspace
wsf= r'F:Mosaic 2000\Reflectance'
rasters = glob.glob(os.path.join(ws,'*.tif'))
for raster in rasters:
try:
    '''bands1 = os.path.basename(raster) + '\Band_1'
    Refl1 = arcpy.sa.Int(arcpy.sa.Int(arcpy.sa.Int((arcpy.sa.Raster(bands1))- 1) /  0.0508))
    outname = os.path.basename(raster)[15:23]+'_B1.tif'
    Refl1.save(os.path.join(wsf,outname))
    print outname
except Exception as e:
    arcpy.AddError(e.message)

try:
    bands2 = os.path.basename(raster) + '\Band_2'
    Refl2 = arcpy.sa.Int(arcpy.sa.Int(arcpy.sa.Int((arcpy.sa.Raster(bands2))- 1) / 0.0254))
    outname = os.path.basename(raster)[15:23]+'_B2.tif'
    Refl2.save(os.path.join(wsf,outname))
    print outname
except Exception as e:
    arcpy.AddError(e.message)'''

try:
    bands3 = os.path.basename(raster) + '\Band_3'
    Refl3 = arcpy.sa.Int(arcpy.sa.Int(arcpy.sa.Int((arcpy.sa.Raster(bands3))- 1) / 0.0363))
    outname = os.path.basename(raster)[15:23]+'_B3.tif'
    Refl3.save(os.path.join(wsf,outname))
    print outname
except Exception as e:
    arcpy.AddError(e.message)

try:
    bands4 = os.path.basename(raster) + '\Band_4'
    Refl4 = arcpy.sa.Int(arcpy.sa.Int(arcpy.sa.Int((arcpy.sa.Raster(bands4))- 1) / 0.0423))
    outname = os.path.basename(raster)[15:23]+'_B4.tif'
    Refl4.save(os.path.join(wsf,outname))
    print outname
except Exception as e:
    arcpy.AddError(e.message)


texts: ArcGIS REST API - return statistics grouped by year

texts:
I would like to query a feature service and retrieve statistics grouped by year. The problem is I don't have a year field, only a date field. Is there a way to convert the date field into a year field? It needs to be live updating since this will be used in an app and a feature service that is constantly getting new records.
My current python code that works with an existing group-by field:
queryParams = {
    'f': 'json',
    'where': "1=1",
    'outStatistics': "[{'statisticType': 'sum','onStatisticField': 'RDM_lbsac', 'outStatisticFieldName': 'RDM'}]",
    'groupByFieldsForStatistics':'Plot_Number'
}
response = requests.post(url = fsURL, params = queryParams)

I could load in a spatial dataframe using the Python API, but the data set is large and this would take a lot of time every time the app loads.


texts:
Possible Duplicate:
  What is the right solution (and map projection) to compute distances between points located all over the World?

I want to know how I could calculate the center-to-center distance between two (or even more) circles, based on their center's Geo-coordinates values?

I'm new in this area, gonna use this one for Goolge Maps API 3 via Javascript and PHP. Probably the algorithm should be enough, however more details are appreciated. :)

Update:

The circles are not necessary to mention as part of the problem, as it could apply to two points, instead of two circles' centers. Actually I'm trying to do some overlap detection stuff on a crisis mapping application.

The locations are relatively close to each other, however I prefer to go for a fail-safe solution which could apply on most of the cases (if not all).

Also it doesn't need an accurate result, as all the visualizations on this web application are somehow approximate by default, so there is no great rush for having a perfect-accurate solution for now, but same as above, I personally prefer to go for the best choice, if it's not gonna add so much difficulties on the solution.


texts: Using Python to compute the distance between coordinates (lat/long) using haversine formula and print results within .csv

texts:
I have a .csv file with 3 columns. The first column is titled Distance and it is empty. The second and third columns are titled Longitude and Latitude which contain coordinates, for example:

Distance,Longitude,Latitude

,-77.60483,40.31117

,-77.58167,40.3045

,-77.52883,40.24917

,-77.508,40.14917

,-77.49617,40.13117

I have wrote a script that computes the distance between the two points using the haversine formula (see below), but I have to manually enter the coordinates, which is not desirable:

import math

lat1 = 40.31117
long1 = -77.60483
lat2 = 40.3045
long2 = -77.58167

def distance_on_unit_sphere(lat1, long1, lat2, long2):

# Converts lat & long to spherical coordinates in radians.
degrees_to_radians = math.pi/180.0

# phi = 90 - latitude
phi1 = (90.0 - lat1)*degrees_to_radians
phi2 = (90.0 - lat2)*degrees_to_radians

# theta = longitude
theta1 = long1*degrees_to_radians
theta2 = long2*degrees_to_radians

# Compute the spherical distance from spherical coordinates.
# For two locations in spherical coordinates:
# (1, theta, phi) and (1, theta', phi')cosine( arc length ) =
# sin phi sin phi' cos(theta-theta') + cos phi cos phi' distance = rho * arc    length

cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) +
   math.cos(phi1)*math.cos(phi2))
arc = math.acos(cos)*6371 #radius of the earth in km

return arc

print distance_on_unit_sphere(lat1,long1,lat2,long2)

My goal is to populate the empty distance column using Python by computing the distance between the 1st location (-77.60483,40.31117) and the 2nd location (-77.58167,40.3045) and then print that value into the Distance column and then compute the distance between the 2nd (-77.58167,40.3045) and 3rd point (-77.52883,40.24917) and then print that value below the previously computed value in Distance column and so on for the subsequent coordinates. I need some help or an example provided to input the .csv file, run the haversine formula, and then write the result into the appropriate location in an output .csv.


texts: Reprojecting and saving shapefile in gdal

texts:
I have a shapefile I want to reproject based on another projection, and then save the results.

I am trying to do this like so:

from osgeo import ogr, osr, gdal

#tif with projections I want
tif = gdal.Open("C:/path_to_tif/file.tif")

#shapefile with the from projection
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource =   driver.Open("C:/path_to_shape/file.shp", 1)
layer = dataSource.GetLayer()

#set spatial reference and transformation
sourceprj = layer.GetSpatialRef()
targetprj = osr.SpatialReference(wkt = tif.GetProjection())
transform = osr.CoordinateTransformation(sourceprj, targetprj)

#apply transformation
for feature in layer:
    transformed = feature.GetGeometryRef()
    transformed.Transform(transformed)

    #save shapefile
    to_fill = ogr.GetDriverByName("Esri Shapefile")
    ds = to_fill.CreateDataSource("C:/path_to_projected/projected.shp")
    outlayer = ds.CreateLayer('', None, ogr.wkbPolygon)
    outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
    defn = layer.GetLayerDefn()
    feat = ogr.Feature(defn)
    feat.SetField('id', 123)

    geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb())
    feat.SetGeometry(geom)
    outlayer.CreateFeature(feat)

but the resulting shapefile is empty and feat.SetGeometry(geom) returns 0


texts: Why empty tiles when zoomed "too much" in Leaflet.js?

texts:
If I set maxZoom to any high number like 100, and then zoom in beyond the maximum "zoom level" for the actual tile layer, Leaflet.js doesn't do the expected and show the highest resolution tile, stretched out, but instead just shows blank tiles.

Can it be changed somehow so that it keeps using the highest resolution and stretches it out for further zooming?

I know the details don't get better, but that's not why I want to keep zooming.


texts: Python-gdal write a GeoTiff with binary color and NaN

texts:
I am generating a simple two-class (binary) geotiff from a numpy array with 3 values: 1, 2, and  NaN. I wanted to display the geotiff with distinct colors so I used the color table. However, color table only supports the Byte or UInt16 datatype, which will convert the NaN to zeros.

Is there a way to write a binary geotiff including the NaN with colors?

filename = 'test.tif'
ds = gdal.Open(filename, gdal.GA_ReadOnly)
outfile = 'classify.tif'
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outfile, cols, rows, 1, gdal.GDT_Byte)
outdata.SetGeoTransform(ds.GetGeoTransform())
outdata.SetProjection(ds.GetProjection())
outband = outdata.GetRasterBand(1)
ct = gdal.ColorTable()
ct.SetColorEntry(1, (0,0,102,255))
ct.SetColorEntry(2, (0,255,255,255))
outband.SetColorTable(ct)
outband.WriteArray(result)
outband.FlushCache()
outdata = None
outband = None
ds = None


texts: esri draw graphic at a fixed screen coordinate

texts:
I am using ESRI silverlight. I have graphic objects on the screen and each graphic object has an attached silverlight label. What I want to do is to have each label shown at a fixed screen (pixel) distance from graphic object. This seems not possible as esri draws the labels according to the envelope given in map coordinates.

Is there an easy way to show labels at fixed screen coordinates, not in map coordinates?


texts: QGIS: Edit attribute values of GeoJSON file: toggle editing greyed out

texts:
I'm trying to edit the values of attributes in QGIS, but it does not work. I tried to google an answer, but couldn't find anything so far. This is what I'm doing:

Open QGIS.
Drag in a GeoJSON file. It becomes a vector layer.
Menu bar > Layer > Toggle Editing.

I cannot toggle editing. It is greyed out. Am I missing something (I'm relatively new to QGIS)?

QGGIS Version

QGIS version
2.10.1-Pisa
QGIS code revision

Compiled against Qt
4.8.6
Running against Qt
4.8.6
Compiled against GDAL/OGR
1.11.2
Running against GDAL/OGR
1.11.2
Compiled against GEOS
3.4.2-CAPI-1.8.2
Running against GEOS
3.4.2-CAPI-1.8.2 r3921
PostgreSQL Client Version
9.4.4
SpatiaLite Version
4.2.0
QWT Version
6.0.2
PROJ.4 Version
491
QScintilla2 Version
2.9


texts: Build Network after changing streets feature class in ModelBuilder?

texts:
Previous posts have said a network cannot be created in ModelBuilder. However, I have been able to create a network after going through part of the wizard, but stop short of build network, then go back and edit the streets network.

For some reason, when I iterate through the creation of new streets, I cannot input those new streets into the model. I get an error 99999 and it says the problem is that my "OBJECTID" is all null values, which is the same in my previous work outside of  MB that worked.

How can I execute "build network" after editing the streets feature class?


texts:
I am using QGIS and am looking to get started on creating a geodatabase. I have point, line and polygon data. I want to arrange my database in categories that can store various point, line and polygon.  I currently just have a big collection of shape files which I know is not the way forwards!
The amount of processing I need to do with this data is limited, it is mostly just for viewing.  Sometimes I need to do basic processing such as buffer, and clip functions. In the future, I would be hoping that this database can also store photo data with my point data.
Back in the day in school I learned ESRI geodatabase.  I have had some issues with the QGIS geopackage which is similar and am hoping to move to something different.

I am looking to move onto either an SQLite or PostGIS system but I'm just not sure where to start.  Once I decide which one to use I have found a lot of online resources where I can teach myself how to get it set up.  At this point, I am  just looking for a recommendation of which system to dive into.

The datasets that I work with are simple, however, they contain large amounts of data.


texts:
I've just installed GDAL, but I don't know how to import/use ogr2ogr.

I've tried:

from osgeo import gdal
from osgeo import ogr

and searched both for ogr2ogr, but it wasn't there.


texts:
I'm trying to get a weather image to display on Google Maps, however in my comparisons to others it seems as if my image is slightly to small (Vertically).. I'm wondering what I might have done wrong in the process..

I have my latitude / Longitude coordinates

$NWlat = '50.406626367301044';
$NWlng = '-127.620375523875420';
$SElat = '21.652538062803';
$SElng = '-66.517937876818';

And when converting them to pixels on Google maps (at zoom level 5) I get:

1191,2763 and 2582,3591

Taking the diferences between these I get:

1391 x 828

I take my image which is currently sized at 3400x1600 and force it into that size above (doing this in photoshop). I then use the Crazed Monkey tile cutter to cut the tiles and place them on the map.. things are slightly not verticle enough.

CMD:

googletilecutter.sh -z 5 -o 5 -t 4,10 -p 167,203

See to two images below, and how the graphics are slightly not as elongated from wrong to right. (take note of location of the blue around kentucky and over by rhode island)

I would greatly appreciate any help with this as I've been struggling now for a bit.

Thank you so much

Wrong:

Right:


texts: QGIS automatic fill columns from referenced csv

texts:
I have what is probably quite an easy problem, however I'm not able to figure it out by myself. I work as a field researcher doing bird inventory. I collect the data via QGIS on an android-tablet, which works fine. For the later evaluation and display of the collected data, I need certain columns in the produced shapefile, which I want to be (at least partly) filled out automatically.
There is one column with the bird common names, another on with the scientific names and a third one with abbreviations.
What I do in the field is fill out the common names by drop down menu (value map) from a csv table, which also contains the scientific names and abbreviations.

The csv looks like the following example:
Common name;Scientiffic name; Abbreviation
Golden Plower;Pluvialis apricaria;Glpl

Later I complete the other two columns at the office.
What I now would like to do is to get the columns, which shall contain the scientific names and abbreviations, to be filled automatically by the associated values from the mentioned csv.
Does anybody knows a way to do that? It would make work way much easier for me and save lots of time during the evaluation of the data...


texts: Adding two to every value of field using ArcMap Field Calculator

texts:
I am working on a road network and each section is assigned a certain value and increases by 1 each section until it gets to around 500 which is the end.

A change I had to make was adding two sections to the beginning of the road which now means every value on this segment is two below what it should be with the addition of these two sections.

I know there's a way to do this in field calculator, but it is definitely my weakest area of expertise in ArcMap.

What kind of function should I use to make these changes?


texts: Automating Select by Location in ArcGIS Desktop?

texts:
I have 2 shapefiles (A and B) each with thousands of polygons. I want to select polygons from A which intersect a specific polygon from B, and save the selected polygons from A as a new shapefile.

I can do this using the 'select by location' feature and 'create layer from selected features'.

However, I want to repeat this for every single polygon in B.

Is it possible to automate such a task? I don't know Python, but may be able to have somebody do this for me if that is the best method.


texts: Loading a GML file to PostGIS using Python and GDAL vectortranslate: encoding problems

texts:
I have some GML files I want to load into my PostGIS database using the GMLAS reader in GDAL in a Python script. The script is run from a batch file. For example this dataset from Geonorge - https://kartkatalog.geonorge.no/metadata/grunnskoler/db4b872f-264d-434c-9574-57232f1e90d2.
I can import a GML file to my PostGIS database using ogr2ogr from the osgeo4w shell, and it goes smoothly, it does not complain about encoding and the Norwegian characters that are in table names and column names, and everything is imported correctly.
ogr2ogr -f PostgreSQL PG:"host=host port=port dbname=dbname user=user password=pass" GMLAS:"<path to gml>" -nlt CONVERT_TO_LINEAR -lco SCHEMA=output_schema

However, when I try to do the same with gdal.VectorTranslate in a Python script, I get problems with encoding. "invalid byte sequence for encoding "UTF8":"
srcDS = gdal.OpenEx('GMLAS:path_to_gml', gdal.OF_VECTOR, )

gdal.VectorTranslate(
        postgisDS,
        srcDS,
        format = 'PostgreSQL',
        geometryType = 'CONVERT_TO_LINEAR',
        layerCreationOptions=['SCHEMA=output_schema']
        )

My GML file says it is in UTF8, the first line is <?xml version="1.0" encoding="UTF-8"?>. Opening the file in Notepad++ and Notepad says the file is in UTF8.
My database has UTF8 encoding. I am doing this on a Windows Server 2012 R2 Datacenter machine using Python 3 and gdal installed from osgeo4w. Gdal version is 3.1.4 og Python is 3.7.
I have read of different posts and websites on this, and I cannot seem to get anything to work. I have tried setting PGCLIENTENCODING, and changing the codepage of the command line, and lots of other stuff. I am now very confused, and still can’t get it to work.
The script works fine on a Linux machine I have access to, but I need it to work on Windows.
Any ideas?


texts: Export geotiff as pdf using external python console

texts: I have a geotiff which I need to export as png or pdf using python only, not the python console inside QGIS. So the whole procedure of opening the PrintComposer, ideally adding colorbar and adjusting colormap, and exporting to pdf (or any image format) needs to be implemented in an external python script. The whole procedure needs to be fully automated, which is why it does not help to simply export from existing PrintComposer, for which there are plenty of questions already, at least explaining how to do it from the QGIS python console.


texts: What surface's area is used to compute area of country?

texts:
Premise:

Earth's surface geometry is modelled by means of different geoids -depending on the objectives. The roughest and most basic are ellipsoids surfaces. Yet, on a small scale, e.g., area of private parcels of farmer's land, municipalities areas - direct measurements, local geometry and neglecting any curvature are used to compute areas. At this scale, 100 m by 100 m parcel does have an area of 10'000 sqm - even if it lays on a slope 45° inclined.

Problem

If we increase the scale and consider a country like Switzerland - in which a good part of the territory is made up of slopes, it is not clear to me how is the country area calculated. In other words:

Is the area obtained by a bottom up approach - summing up the small areas of the small parcels of territory? in which case the overall area accounts for the effective, real, extension of the land …
Is the area calculated by calculating the area of a simplified country surface which is NOT the union of all the small parcels? in which case the country area is NOT equal to the sum of the areas of the small parts …
If 2 is the case, does the simplified country surface still include slopes or is it made "flat" (slope's area is the area of the projection on the base)?
ACID Test: if a country includes great many low inclined slopes would its area change if the slopes were flattened?


texts: Label orientation in Mapbox Studio or Mapbox GL JS

texts:
Is it possible to make a label orientation a function of the shape of the polygon it is labeling?

I imagine it's possible with Mapbox GL JS with https://www.mapbox.com/mapbox-gl-style-spec/#layout-symbol-text-rotate, but that seems like it might be a non-trivial JavaScript task (?). The end result will be a Mapbox GL served vector basemap.

Medicine Lake default orientation:

Medicine Lake ideal orientation (but I don't necessarily want any other lake label rotated):


texts:
I get the following error when importing GeoPandas. I using pip and anaconda to install GeoPandas, but I still get the same error. What I should do?
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-9-b86cb47439ab> in <module>
      4 import numpy as np
      5 from shapely.geometry import mapping
----> 6 import geopandas as gpd
      7 import rasterio as rio
      8 from rasterio.plot import plotting_extent

~/opt/anaconda3/lib/python3.7/site-packages/geopandas/__init__.py in <module>
      1 from geopandas._config import options  # noqa
      2
----> 3 from geopandas.geoseries import GeoSeries  # noqa
      4 from geopandas.geodataframe import GeoDataFrame  # noqa
      5 from geopandas.array import points_from_xy  # noqa

~/opt/anaconda3/lib/python3.7/site-packages/geopandas/geoseries.py in <module>
      7 from pandas.core.internals import SingleBlockManager
      8
----> 9 from pyproj import CRS, Transformer
     10 from shapely.geometry.base import BaseGeometry
     11

ImportError: cannot import name 'CRS' from 'pyproj' (/Users/rohullahnajibi/opt/anaconda3/lib/python3.7/site-packages/pyproj/__init__.py)


texts: PostgreSQL/PostGIS error "pseudo-type record"

texts:
I've got, in PostgreSQL/PostGIS,

create table datos (
gbifID int primary key,
orden varchar (50),
eventDate varchar (50),
geom geometry);

When I run

select distinct (geom,eventdate), geom
from datos
where orden='x';

it works OK. However,

create view esfuerzo as
select distinct (geom,eventdate), geom
from datos
where orden='x';

returns
ERROR:  column "row" has pseudo-type record
Estado SQL: 42P16. How can I create such a view?


texts: Calculating fragmentation parameters at high resolution

texts:
I'm trying to evaluate fragmentation of an ecosystem. While It is easy to calculate in ArcGIS 10.1 values such as Area, perimeter or minimum distance between fragments, I'm having problems calculating other parameters that are used in software such as Fragstats. I tried to use this software but I have problems with the memory, as my raster have a very high resolution.

I wanted to calculate first the "edge effect" also called Core Area in Fragstats. I mean, if I subtract for example the first 100 m from the border of those fragments, how much area do remain? I tried to make this with buffers but it is very difficult to calculate polygon by polygon.

Is there any way to calculate this more automatically? Maybe another options to calculate similar parameters to those of Fragstats in ArcGIS?


texts: Converting GDAL dataset into numpy array

texts:
I am trying to do data augmentation on sentinel 2. Everything else is working the only problem I have is the output I am getting is only one band when the input image is with 10 band. I also tried with rasterio but the result are still only one band. Can someone tell me why i am only getting one band as ouptut. Thanks
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
# import io
from PIL import Image
from osgeo import gdal
import cv2
import rasterio
import numpy as np
from osgeo import gdal

datagen = ImageDataGenerator(
        rotation_range=180,
        width_shift_range=0,
        height_shift_range=0,
        shear_range=0,
        zoom_range=0.5,

        horizontal_flip=True,
        vertical_flip = True,
        fill_mode='nearest')    #Also try nearest, constant, reflect, wrap

a = 1
while True:
    ds = gdal.Open("E:\\opencv\\12.tif")
    myarray = np.array(ds.ReadAsArray())
       ""Rasterio code""
    #fp = r'E:\\opencv\\12.tif'
    #img = rasterio.open(fp)
    #img = img.read()
    """"""
    x = img_to_array(myarray)
    x = x.reshape(x.shape + (1, ))
    i = 1
    for batch in datagen.flow(x, batch_size=16,
        save_to_dir='E:\\opencv\\',
        save_prefix= a,
        save_format='tif'):
        i += 1
        if i > 1:
            break
    a += 1
    if a > 3:
        break


texts: Can't find the toolbar for sorting the table of contents in alphabetical order in Arcmap!

texts: I am using ArcGIS Desctop 10.0 Service Pack3 and I'd like to sort the table of contents in alphabetical order but following Sorting the table of contents in alphabetical order I can't find the "Production TOC Manager" toolbar. Am I missing anything or do I need to activate anything?


texts:
I am using ArcGIS Desktop to view one shapefile.

Only a select number of features are being displayed. The same number are showing in the attribute table even though the database file (.dbf) has a complete list with many more features.

I then tried opening the shapefile in QGIS and it informs me that the file is invalid. Upon inspection of the shapefile, I found that the .shx file is missing.

Can I recreate the missing .shx file so that all the attributes/features can be loaded?


texts: Estimate forest stand height from Landsat imagery?

texts:
Is it possible to estimate the height of forest stands using optical (preferably Landsat) imagery?

I came across this article while researching this, but it's over 20 years old and there doesn't appear to be much in the way of more recent literature. Can spectral responses be used for deriving forest structural characteristics, or is this asking too much of the data?

There is some information on combining LIDAR data and optical imagery, which makes sense, but I'm looking for an approach that can be applied at a broader scale without a huge cost.


texts: Turning vector grid into individual shapefiles using QGIS?

texts:
I am using QGIS and I created a vector grid on a crop field in order to cut it into plots. I want to save each grid cell into individual polygon shapefiles automatically.

Here is an example of a vector grid I would like to split into 117 shapefiles:

Eventually, I did it manually with the selection tool and "Save as" but it is really time-consuming when there are more than 100 plots to save.

Is there any way to do it faster and more automatically?

Each plot as an attribute named "Plot" that can discriminate against them (basically being id + 1). It's the attribute used here to label them.


texts: Generating transects perpendicular to track generated from a PPK GPX file

texts: I am new to QGIS and barely getting the hang of it so this question is quite simple. I have a route in a GPX file a generated from a PPK I took along the shore during a high tide. I would like to generate transects perpendicular to this track with a spacing of, lets say every 5 m and a length of 20 m each. Any suggestions as to what pluggin or tool I may use to delimit this area?


texts: Merging attributes of point features at same location

texts:
I have a CSV-file with all the facilities in an area, some sites have multiple facilities so at some locations there will be multiple points e.g. sports hall and swimming pool.
Is there a way to join all the points at one location (they have a unique ID) but retain the information so that if I click on the point it will show that there is a sports hall and a swimming pool not two points one for sports hall and one for swimming pool?


texts: Setting ArcGIS field access restrictions?

texts:
I am using ArcGIS 10.5 and creating features class in a PostgreSQL 9.5 database. But some feature classes have specific fields that some users should not edit. For example I have a "gas pipes" feature class and it has a facility_id column. Admin users can edit this field but Viewer users can not edit it from ArcMap. But I could not restrict the field of feature class. Or I can not hide facility_id field.

Is this possible?


texts: How to get BBOX from WMS GetCapabilities response?

texts:
I have a WMS layer and want use zoomToExtent after CQL filter. So I get the layer BBOX from the GeoServer WMS using a GetCapabilities request.

wms = new OpenLayers.Format.WMSCapabilities();
OpenLayers.Request.GET({
  url:"http://localhost:8080/geoserver/wms/filedata?request=GetCapabilities",
  success: function(e){
    var response = wms.read(e.responseText);
    var capability = response.capability;
    for (var i=0, len=capability.layers.length; i<len; i+=1) {
      var layerObj = capability.layers[i];
      if (layerObj.name === 'filedata') {
        map.zoomToExtent(OpenLayers.Bounds.fromArray(layerObj.llbbox));
        break;
      }
    }
  }
});

But zoomToExtent doesn't work with lon lat, so how to get the BBOX in meters?

My map showed in epsg900913.


texts: Do I need an Add-in or a Stand-Alone Application for a custom form?

texts:
The type of application I want to build for ArcMap is to be run from inside of ArcMap and on a click of a button a custom form appears in front of User in which User enters certain parameters and then by some processing, the application returns resulting image to user in the same ArcMap Document...
So will this application be an Add-in or Stand-Alone Application..??


texts: How to remove a vertex from a MultiPolygon and return the bounding polygon in JTS

texts:
I'm trying to complete a code implementing the Kirkpatrick algorithm in Java using JTS and it's literally the first time doing working with such tools. There's a part that I need to remove a known point P (in the independent set) from a region which is consisted of multiple polygons all intersecting in P and return the remaining hole to be triangulated. The question is how to do remove a point P from a MultiPolygon?
 /**
 * Find a bounding polygon if point p is removed from affected regions
 *
 * @param p
 * @param affectedRegions
 * @return
 */
public Polygon findBoundingPolygon(Coordinate p, List<Geometry> affectedRegions) {
    // TODO: write a function that finds a bounding polygon
    GeometryFactory fact = new GeometryFactory();
    Geometry regions = fact.buildGeometry(affectedRegions);

    // TODO: bound <- a concave hull of regions

    LinearRing shell = fact.createLinearRing(coordsToLinearRingCoords(bound));
    return new Polygon(shell, null, fact);
}


texts: Slope Classification for HRU in a watershed (SWAT Model)

texts:
For the definition of HRU (Hydrological response units), the parameter slope can be classified into maximum of 5 classes. Depending on the definition the results will vary.

Is there a commonly used method to decide on the intervals of such classifications?

In literature, many different classifications are used. Most of the times, no information on how the classes were created is given.

How can I find out what makes sense the most?

Example FAO

0-7
7-15
15-20
20-25
25-30
30-100

Example ArcGIS

In ArcGIS, the slope can be classified with different algorithms. For my watershed, the Natural Invertall (Jenks-Caspell-Algoritm) creates the classes

0-5
5-12
12-20
20-30
30-100

Standard deviation

0-4
4-13
13-22
22-31
31-100


texts: Using LiDAR to identify coniferous and deciduous trees

texts: Is it possible to use only LiDAR data to distinguish between mostly coniferous and deciduous forest? I have a LiDAR data from spring and summer, but they have different point spacing.


texts: Get all points with in some range of a point (or range of point) in MySQL

texts:
I have got a list of points and a fixed range (10m). So what I have to do is for each point I have to find all the points present in a table with the 10m range of the given point.

For example:
Let say I have following list of points:
P1, P2, P3, P4, P5, P6

Now I have to loop through each point.
First I'll get P1 and I need all points (in a table) present inside the circle where radius is 10m and center is P1.
Similarly I'll do for all other points from P2 to P6. My result will be all the outcome points in each loop.

I am not sure whether I have explained well or not. But I can clarify if you can comment.

I also need to improve the solution; where I can calculate the outcome points one at time (by creating a polygon maybe) and then somehow consider the range.

Thanks


texts: Combine attributes (lines) within an attribute table QGIS

texts:
I have an attribute table with 100 LINES.
The table has a DATE field.
The 100 LINES have several DATE values.
I would like to combine all LINES with the same DATE value.

How can i do that?


texts:
pgRoutingLayer does not show the database connected.

QGIS: 2.14.15 and pgAdmin III

CREATE EXTENSION postgis; (ok)

CREATE EXTENSION pgrouting; (ok)

SELECT pgr_version(); (2.0.0,pgrouting-2.0.0,0,d6ed2cb,master,1.53.0) (ok)

How to solve this?


texts: how to use mapproxy to serve WMS from reprojected google maps tiles, i have searched for docs, but there's no docs which describe how to add google maps as source.


texts:
I'm working with the last version of QGIS (3.8.3-Zanzibar), and I'm trying to represent a Rule-Based Symbology from a CSV file with the following rules:

And the data from the prcptot column is this:

As you can see, all the prcptot values are negative (smaller than 0), but the output that QGIS is giving to me is the next:

All the red triangles should be inverted blue triangles... I think that QGIS is taking all the values as positive numbers, but I don't know why or what to do.


texts:
my plugin creates an ASC file that I can load in QGIS as raster.

Sometimes I need to load the raster with an additional rotation.

The only solution I can think is to warp the raster by a custom CRS with an additional rotation.

Is there another smartest way?

EDIT
This is the numpy array that I have. Actually the asc file comes from this array, so all the transformations will be performed on this object.

The code below does not work with this dataset:

import numpy as np
data = np.array([[ 40.,  40.,  40.,  40.,  40.,  40.,  40.],
[ 36.,  36.,  36.,  36.,  36.,  36.,  36.],
[ 32.,  32.,  32.,  32.,  32.,  32.,  32.],
[ 28.,  28.,  28.,  28.,  28.,  28.,  28.],
[ 24.,  24.,  24.,  24.,  24.,  24.,  24.],
[ 20.,  20.,  20.,  20.,  20.,  20.,  20.]])


texts: Using ArcGIS Desktop to calculate Conversion Probability certain land use to other land use?

texts:
I am doing dissertation about [predicting future urban expansion for Shannon River Basin In Ireland using Cellular Automata model and arcGIS] so at this stage I would like to calculation the conversion probability factor (Land use transitional probabilities), as it is advised that the value of this probability can be generated fixed by overlaying two land use maps, so it can see the value of conversion rate of a certain land use to other land use. so that I enter it into the model.

I am kinda stuck on which arcGIS tool of OVERLAYING should I use (union, erase,identify etc), knowing that I have two landuse/landcover maps from EPA website ready with almost 32 different landuse classes for the years 2000 and 2006
and how can I actually calculate the conversion rate ?


texts:
I am novice mapper. I ran into some trouble with QGIS today when trying to add additional census tracts to a layer filter.

The tract exists within the layer I'm filtering, but there seems to be a hard cap of 30 for number of times I can use the OR statement in the filter. I've rearranged the tracts and flipped which ones I am using, but can never get past 30 selections. For reference, my filter query that I ended up with is below: to replicate, use the Census Bureau tract map of Michigan.

"COUNTYFP" = '077'  AND ("NAME" =  '1' OR "NAME" = '10' OR "NAME" = '11' OR "NAME" = '12' OR "NAME" = '2.01' OR "NAME" = '2.02' OR "NAME" = '3' OR "NAME" = '5' OR "NAME" = '6' OR "NAME" = '9' OR "NAME" = '15.01' OR "NAME" = '15.02' OR "NAME" = '15.03' OR "NAME" = '15.04' OR "NAME" = '15.06' OR "NAME" = '15.07' OR "NAME" = '16.01' OR "NAME" = '16.03' OR "NAME" = '16.04' OR "NAME" = '17.01' OR "NAME" = '17.02' OR "NAME" = '18.01' OR "NAME" = '18.02' OR "NAME" = '18.03' OR "NAME" = '29.03' OR "NAME" = '29.04' OR "NAME" = '29.05' OR "NAME" = '30.04' OR "NAME" = '55.01' OR "NAME" = '55.02')


texts: Finding unconnected polylines and connecting them using QGIS?

texts:
I have a pipeline network, during manual digitization have a lot pipes unconnect.

Is there any way to find pipes unconnect and edit it in QGIS?


texts: Command Line alternative to using QGIS GUI to save a Raster as a Rendered Image

texts:
I have a raster image in PNG format which I can open in QGIS but not in Windows Picture Viewer. To do that i figured out that I have to save my raster as a 'Rendered Image' from the 'Save As' (option by right-clicking on it in QGIS GUI).

So great then I could now open my rendered PNG in Windows Picture Viewer.

"Now I need a little help with as to how to do this (i.e. save the Raster as Rendered Image) without the QGIS GUI and via a Command Line"


texts: Import info from mdb file into a project?

texts:
I need to import info on a QGIS project.

The info are in a MDB file (population, mean age, number of cars and so on: one record per region) that is written to be connected with a region shape from an unique ID (the ID is the same in the MDB and in the QGIS shapefile).


texts: How do I generate a correlogram from a GeoTiff using R?

texts:
I'm trying to evaluate the scale of spatial autocorrelation for a bunch of rasters I have as GeoTIFFs.  I'd like to be able to do this entirely using R, in order to maintain smooth workflow.

What's the best way of going from GeoTIFF to correlogram using R?


texts:
I try to make autosaving on adding feature and use next code:

def autosave():
    iface.activeLayer().commitChanges()
    iface.activeLayer().startEditing()
    iface.actionAddFeature().trigger()
iface.activeLayer().featureAdded.connect(autosave)

But commitChanges() execute a featureAdded and I get a loop. Also I try to use attribute tables validation on search not saved features, like this:

def autosave():
    features = iface.activeLayer().getFeatures()
    gid = next(features)[0]
    if "nextval" in str(gid):
        iface.activeLayer().commitChanges()
        iface.activeLayer().startEditing()
        iface.actionAddFeature().trigger()

iface.activeLayer().featureAdded.connect(autosave)

And this code executes loop too. How can I avoid this problem?


texts: How to select polygons which are situated between two points using select by location?

texts: I am working with ArcMap 10.3. and have a polygon layer and a point layer. I've selected two point layers and now I want to select the polygons which are situated between those two points. Is there a chance how I can do that with "select by location"?


texts: Enable both function on click and highlight on hover in Leaflet+cartoDB

texts:
Currently, on http://velo300.free.fr/debug,

feature highlights on 'featureOver' creating a white new feature on top of mouse-overed feature.

Problem: i need also to work with 'featureclick' for another function.

How can i manage to create this new white feature below the mouse-overed feature and not on top?

Other relevent solution: how to use setcartoCSS() wihtout redraw rhe whole layer?


texts: How do I delete duplicates based on multiple fields

texts: I want to delete duplicates based on my time/date column (1st column) but I want to keep the ones that have a 1 in the 2nd column. Essentially, I want to delete all rows that are a duplicate of 1st column AND are "2" in the 2nd column. How do I do this in ArcGIS?


texts: Generating equidistant Points inside box which also contains MultiPolygon

texts:
I have a box bounded by min_x, min_y, max_x, max_y, I also have a multipolygon = poly contained in that box. I'm trying to equidistantly scatter points throughout that box, outside of that multipolygon, for the eventual purpose of using a modified A* algorithm to find a route between two points.
Right now, I can generate points as such:
from shapely.geometry import Point
import numpy as np

x_range = np.linspace(min_y,max_y,round((max_y - min_y)*2.5),False)
y_range = np.linspace(min_x,max_x,round((max_x - min_x)*2.5),False)

points = []
point_dict = dict()
i = 0

for x in x_range:
    j = 0
    for y in y_range:
        p = Point(y,x)
        if not poly.contains(p):
            points.append(p)
            point_dict[i,j] = p
        j+=1
    i+=1

They're in an x-y grid along the map, but I think a more efficient/natural way to generate points for a path would be to create this grid then add a repulsion factor to each point so they're all pushed away from each other inversely proportional to their distance AND from the multipolygon and the bounding box. Eventually I imagine this would end up in some state equilibrium.
I stumbled upon the following question: Electric force between particles using numpy arrays.
And was able to get an idea for how someone might approach this problem, but I'm still unsure how I could bound points to keep themm inside of a box and outside of the polygon/multipolygon that the box also contains.


texts: I install QGIS 2.99 development version on my windows 10. I see in the installation folder I have folders python27 and python36. I want to ask can i use python 3.6 to take QGIS python API in standalone app? Does someone have experience withpython 3.6andQGIS 2.99` to help me? Are there some steps with others version QGIS to take python api ?


texts: Regular expression to extract road labels from OSM other_tags ref

texts:
Extended question to the issue: Regular expression substring for labelling

The answers I've found on the link above were fine for "simple" road labels.

Unfortunately, I don't have a big SQL-knowledge, therefore, I would like to extend Erik's question.
How can we extract other road labels (alphanumeric) from the OSM lines layer "other_tags" between the ("ref"=>") and the last (")?

Examples:

"ref"=>"AP 520" -> AP 520
"ref"=>"M123-4" -> M123-4
"ref"=>"EN 13-5" -> EN 13-5
"ref"=>"AG-52" -> AG-52

In order to avoid to write thousands of different expressions depending on each example, isn't any expression that extracts the content only of "ref"=>"xxxxxxxx"?

I tried everything using the "regexp_substr".


texts: Why does polygon.intersection(complex_linestring) split the linestring and mess up Z?

texts:
I have a complex (self crossing) line with Z values which I'd like to truncate to the area inside a box. The problem is Shapely splits the line into simple (non-crossing) sections and the Z value at the intersection of the loop is incorrect.

My high-level goal is taking a pile of GPS tracks, trimming them to only include the area in a box and calculate speed within the box. I'm using Z to represent seconds. I suspect this problem is deep in geos. Can you suggest a better approach?

I made a pretty picture of this example at http://i.stack.imgur.com/pQxtH.png but it may not work since I'm new here.

from shapely.geometry import LineString
import shapely.geometry
box = shapely.geometry.box(1, 1, 4, 5)
line = LineString([(2, 0, 0), (2, 4, 4), (3, 4, 5), (3, 3, 6), (0, 3, 9)])
[[p for p in l.coords] for l in line.intersection(box)]

results in:

[[(2.0, 1.0, 1.0), (2.0, 3.0, 5.0)],
 [(2.0, 3.0, 5.0), (2.0, 4.0, 4.0), (3.0, 4.0, 5.0), (3.0, 3.0, 6.0), (2.0, 3.0, 5.0)],
 [(2.0, 3.0, 5.0), (1.0, 3.0, 8.0)]]

but I was expecting only this segment:

[(2,1,1), (2,4,4), (3,4,5), (3,3,6), (1,3,8)]

shapely.__version__
'1.2.17'
shapely.geos.geos_version_string
'3.2.2-CAPI-1.6.2'


texts: Unknown set of coordinates Northing Westing

texts:
I've been given a list of sites of vegetation with a "Northing/Westing" location e.g

52 43.569 - 001 44.833

I want to convert them to Easting Northing so I can load them into ArcMap. I've tried to convert the original coordinates but having a hard time determining the best way.

The best method so far is to convert them to Long/Lat

52.43569N / 1.44833W

in gridreferencefinder. However the location assigned to each record is different to the location found when using these new coordinates, off by about 20-30 miles.

Anyone know how to convert these locations to something more workable?


texts: Converting to EPSG:3857 produces results in the x 100 000s

texts:
Using Google Maps I get coordinates around the hundreds values (in an Australian location). For example:

However when I transform coordinatess from epsg:28356 to epsg:3857 (Google Map's projection), I get results in the ten's of millions.

Does anyone know why this is?


texts: How can I add custom label to overpass turbo result

texts:
For example given this overpass turbo query: http://overpass-turbo.eu/s/ctg

How could I change all of the names shown over those areas to say "foo" instead?

What I have tried that does not work:

Change the text:name to text:foo
Use the set tag=value syntax like this:

{{style:
node,way,relation {
set bar=foo;
text:bar;
}
}}
Add " or ' around the tag and text values


texts:
How to run this plugin from PyQGIS??

http://docs.qgis.org/2.2/en/docs/user_manual/plugins/plugins_interpolation.html


texts: Getting the width or maximum distance of multilinestrings

texts:
I have several multiline strings. Each string consists of 2 to 8 almost parallel lines. I now want to have the largest distance between the two outer lines or calculate the maximum width of my multiline strings.

So I'm looking for something like this:
ST_Width(ST_Boundary(ST_GeomFromEWKT(geometry)))

I have a few question about this:

How are length and width of a boundary defined?
Is this dependent on the position in the room?
Is the larger value automatically the length and the smaller value the width?

In my case all lines are about the same length. However, the length of a line can be shorter than the distance between the innermost and outermost line (case 2).
The orientation in my room can also be different (case 3).

Therefore the following sketch illustrates my problem (in red the value I want to calculate):

I'm overthinking my strategy. The lines of the multi line strings are ordered (ID 0 to 8). For each set of multilines I try to use st_distance between the line with lowest and highest ID. Hopefully, it will work.

Otherwise, I try ST_ConvexHull. It seems a good way to solve my problem, too.


texts:
I have one ESRI Shapefile as a map with many objects in it.

I have another ESRI Shapefile only with borders in it. This shapefile is bigger in size (covering a bigger area) with many border lines. The border lines form many mosaic pieces (tiles).

How can I cut a part out of the map shapefile with the border shapefile as one mosaic piece?

Where is the answer described in the Manual (en or dt) or in the answer to another question?


texts: PostGIS and GeoServer - alternative solution to serve combined spatial and time tables with high resolution data

texts:
This is the follow up of this other question of mine, although you don't really need to see it to understand the problem.
I am trying to figure out what would be the best (maybe most usual) way of serving large tables that have both a spatial and a time dimension.
First of all: what I have
I have two tables (not the actual that I have, but you'll get the point):

precipitaion: store precipitation data at daily (freq = 'daily') and hourly (freq = 'hourly') frequencies for 7 months so far (I receive regular updates).

Column  |         Type
--------+----------------------
 id     | varchar
 date   | date
 start  | time
 end    | time
 freq   | varchar
 prec   | float(4)
(2,481,069 rows)

areas: store multipolygon areas for all Europe down to municipality resolution.

Column  |         Type
--------+----------------------
 id     | varchar
 geom   | geometry
(162,573 rows)

What I need
I am developing an app that lets the user select a day, or a time in one day and shows in a (OpenLayers) map the areas with the amount of precipitation for that selected time.
I am using GeoServer to serve the data with CQL_FILTERs to pick the selected time (see other sections to understand how the data is built).
What I was doing
I created one table view for daily frequency precipitation and one for hourly frequency precipitation, like so:
create or replace view daily_prec
as select * from
    (select
        id,
        "date",
        "start",
        "end",
        sum(prec) as prec,
    from precipitaion
    where freq = 'daily'
    group by "date", "start", "end", id) as prec,
    areas.geom
where prec.id = areas.id;

For the daily_prec I have 4,634,430 rows, and much more for the hourly_prec.
Serving these two views as WMS layers through GeoServer was slow, as all getMap and other requests took more than 1 sec to arrive, so...
...what I am doing now
I created materialized views so that I could create index and spatial index on the tables.
Now, the requests arrive much faster, and I am happy with the overall performance.
The problems are:

they require a lot of disk space (a single materialized view table can take up 60GB)
they are very slow to be created (they take about 40 min).

Is there a better approach?
Given that my data needs frequent updates (also the historical data might change as the way it is calculated might vary because we are in a testing phase), and preparing tiles probably is not the way to go (right now at least), would you recommend a different approach for such situation?
Would you recommend using a different type of database and/or service?
I am using PosgreSQL 10 and PostGIS 2.5.


texts: TIN when generated goes off the limits ARCGIS 9.3

texts:
i have to creat TINs for a specific areas so the original points or countours are exactly whithin the geographical limits of the area, but when i generate TIN it goes out of those limits and its normal based on the principle of triangulation but is there anyway to gard the same shape and limits after creating a TIN ? right now im using "paint" to clean the mess so please help me :)
thanks


texts: Updating table setting a CASE with ST_Intersection

texts:
I add a new geometry column to the "mad" table so I could update it with the polygons from an Intersection with the polygons of the table "clip".

Since my database is huge an it was taking a lot of time to comput I opt for this solution "Getting intersections the faster way" (http://postgis.net/2014/03/14/tip_intersection_faster/). These SELECT worked for me but I wanted to do an UPDATE. I tried using this command:

UPDATE mad SET geom_p=CASE WHEN ST_CoveredBy(mad.geom,clip.geom)
    THEN mad.geom
    ELSE ST_Multi(ST_Intersection(mad.geom,clip.geom)) END AS geom FROM mad AS mad
  INNER JOIN clip AS clip
     ON (ST_Intersects(mad.geom,clip.geom) AND NOT ST_Touches(mad.geom,clip.geom)));

ERROR:  syntax error at or near "AS"
LINE 1: ...(ST_Intersection(mad.geom,b.geom)) END AS geom FR...


texts: I'm working with qgis v. 2.0.1. I have some data points and I just wanted to create a new shapefile with line geometry. In the new vectorlayer-window, I added some attributes (approx. 10) and klicked ok but it didn't create a shapefile (nothing happened). I tried again without adding any attributes and it worked. Why is there a problem when I want to add some attributes? Of course I could add attributes with the Table Manager or manually in the attribute table but I'm just wondering why it's not working this way. Thanks.


texts: MapServer | Cant Get WMS GetFeatureInfo To Return GeoJSON

texts:
I have MapServer serving up a whole bunch of data, and everything is working good.  One area I am stuck is getting the WMS GetFeatureInfo request to return GeoJSON.  Regardless of what I try, I just get plain text back.

The applicable pieces of my mapfile are below.  My understanding is that the OGR/GEOJSON driver (in OUTPUTFORMAT) makes it so that I dont need to use a template.  Is that correct?  I've found that if I dont include some dummy TEMPLATE in the LAYER - the layer is not queryable.  With the below configuration - regardless of what I set info_format (or outputformat) to - I just get plain text response back.

Example response:

  GetFeatureInfo results:

  Layer 'treecrowns'   Feature 42342:

OUTPUTFORMAT
        NAME "geojson"
        DRIVER "OGR/GEOJSON"
        MIMETYPE "application/json; subtype=geojson"
        FORMATOPTION "STORAGE=stream"
        FORMATOPTION "FORM=SIMPLE"
END

 WEB
        QUERYFORMAT "geojson"
        METADATA
            "wms_title" "DAM WMS"
            "wms_abstract" "DAM WMS"
            #"wms_onlineresource" "http://mapserver//mapserv"
            "wms_srs" "EPSG:4326"
            "wms_encoding" "UTF-8"
            "wms_enable_request" "*"
            "wfs_enable_request" "!*"
            "wfs_encoding" "UTF-8"
            "wfs_getfeature_formatlist" "gml,geojson"
        END
END

LAYER
    NAME "treecrowns"
    TYPE POLYGON
    TEMPLATE fooOnlyForWMSGetFeatureInfo # For GetFeatureInfo
    EXTENT -180 -90 180 90
    CONNECTIONTYPE postgis
    PROCESSING "CLOSE_CONNECTION=DEFER" # For performance
    CONNECTION "host=x.x.x.x port=5432 dbname=x user=x password=x"
    DATA "geom from treecrowns"
    METADATA
        "wms_title" "treecrowns" # For WMS
        "wms_srs" "EPSG:4326 EPSG:3857 EPSG:4978" # For WMS

        "wfs_enable_request" "*" # Enable WFS for this layer
        "gml_include_items" "all" # For GetFeatureInfo and WFS GetFeature (QueryBuilder)
        "ows_geom_type" "Polygon" # For returning geometries in GetFeatureInfo
        "ows_geometries" "geom" # For returning geometries in GetFeatureInfo
        "wfs_featureid" "name"

        "wms_metadataurl_href" "http://www.example.com/bar" # For metadata URL
        "wms_metadataurl_format" "text/html" # For metadata URL
        "wms_metadataurl_type" "TC211" # For metadata URL
        "wfs_getfeature_formatlist" "gml,geojson"
    END
    STATUS ON
    PROJECTION
    "init=epsg:4326"
    END
    CLASS
        NAME "treecrowns"
        STYLE
            OUTLINECOLOR 0 0 255
            WIDTH 2
        END
        STYLE
            COLOR 0 0 255
            OPACITY 40
        END
    END
END


texts: How to compress enterprise geodatabase with stored procedure?

texts: I'm looking for a stored procedure in ArcSDE for SQL Server which performs geodatabase compression to include it in our daily maintenance plan. Is there any such thing in ArcSDE 10.0?


texts:
In QGIS (2.14.3.) how do I break polyline features in their consisting (single)line features? I have tried the 'Explode lines' geoprocessing tool, but this destroys the whole polyline into very atomic lines consisting of only start- and endpoint.

This visualizes what I have in mind (offset in start/endpoints only for visualisatin purposes):


texts: How to save 'No data value' to a QGIS layer style properties?

texts:
I have many raster tif files (100+) I found this link to the MultiQml plugin very helpful, from this post: Can I change style properties of several layers simultaneously?

However the 'No data' values does not seem to be saved with the style? Is there a way to do this? Any suggestions are very welcome.


texts: Seeking neighborhood boundaries globally?

texts:
Besides Zillow's creative-commons licensed neighborhood boundaries shapefiles, what sources of neighborhood boundaries are there? Not just in the US, but all over the world.

Free/open source is ideal, but paid sources could be interesting too.


texts: Is it possible to access the Catalog window programmatically in ArcMap?

texts:
Specifically, I would like to get the selected object, similar to the way you do in ArcCatalog through IGxDialog.SelectedObject.

Is there a way to do this?

So far I can get the Catalog DockableWindow but am stumped from there.


texts:
I have 2 shapefiles and after identifying their intersections I'd like to know the ratio of intersecting areas, e.g.

  a given shape S in file1 intersects with shapes (x1,...,xn) in file2 with ?proportions (y1,...,yn) (y_i of shape x_i lies in S).

I do not think that planar areas (e.g. shapely) are appropriate for good estimations. How (what projections) should I use to get the closest approximations to the actual area ratios? What libraries in Python or command line tools allow for calculating areas on hyperplanes (local Earth surface)?


texts: JPEG base map doesn't match shapefile locations

texts: Relative Newbie Question: I have Civil War earthworks tied to the State Plane Coordinate system with a jpeg of the terrain as a 100k scale USGS Quad of the area that is not georeferenced. When I work with the features in CAD, they all have georeferences. But when I export the shapefiles out of CAD and into QGIS, all of the shapefiles come into QGIS where they should be. BUT, the jpeg is not only not under them, but it is 10x larger than it should be and displaced by some absurd difference. Can anyone shed some light on what I'm not doing in order to get them all into one area and at the right size?


texts: Get XY value from TiledMapService Layer Tile

texts:
I have an ESRI TiledMapService layer and I have attached an event handler to the TileLoaded.  Given the level, row, and column, how can I get the envelop of the tile in map coordinates? I'm using the ESRI Silverlight API...

What I am trying to do:
 I have a polygon (in the same coordinates of the TiledMapServiceLayer) and I want to examine each tile as it comes in and cut out the pixels of the tile that intersect with the polygon.

Any help would be great!


texts: I´m trying to do a Least Cost Path with Grass and when I try to use r.drain I can't. It say to me error raster map not found.


texts: Splitting population data proportionally using Identity tool of ArcGIS Desktop

texts:
I am preparing a term project and I have the vector data of the roads and zones layer for to prepare rural access index in Istanbul. My zones layer also has the attribute nufus17 which means population in the 2017 year.

I am stuck in the phase of dividing this attribute proportionally to the areas within a grid that I created via using make feature layer and identity to the zones layer. Because when I used the Identity tool areas perfectly split yet population data are not. It is just copied for each split area.

For example I had a polygon area with 1000 population now I have 2 or more polygon areas with 1000 population.


texts:
I'd like to write an app which uses google maps (via javascript library like gmaps.js).
I'll GET the GPS positions from my server and I'll use google maps+js library to draw them on map.
I don't need to use any api-keys.

Is it free and unlimited or even for commercial usage?


texts: Is there a limitation in quantity of joined tables using spatialite views?

texts:
I've one point layer in qgis and several non spatial tables in a spatialite database. Now I want to create a view as you can see in the following:

CREATE VIEW testansicht10
AS SELECT baumkataster.baum_ID AS ROWID,
baumkataster.geom AS geom,
gattung.gattung AS gattung,
art.art AS art,
sorte.sorte AS sorte,
baumkataster.b_hoehe AS b_hoehe,
baumkataster.kro_durchm AS kro_durchm,
baumkataster.kro_ansatz AS kro_ansatz,
baumkataster.age AS age
FROM baumkataster JOIN gattung ON baumkataster.gattung = gattung.gattung_ID
JOIN art ON baumkataster.art = art.art_ID
JOIN sorte ON baumkataster.sorte=sorte.sorte_ID;

When I do this the result is an empty view table though there are entries in my spatial table. If I exclude one of my joined tables (for instance "sorte") in the view I get the desired result - the view table with the existing entries. Is there a limitation in quantity of joined tables? The view should consist of about 65 coulomns (36 are from joined tables).


texts: ERROR unsupported operands type(s) for *:'int' and 'Results' in for cycle

texts:
I'm trying to do a for cycle to calculate standard deviation and mean for each raster in a list and after that to do a mathematical operation between them. Can someone help me to understand what is the problem?

ws2 = r"C:\Users\Elisa\Model_Builder\passo_5_6_7_8\output\raster_tagliati_su_flowacc"
v_name_sc_img = r"C:\Users\Elisa\Model_Builder\passo_5_6_7_8\output\outliers\{0}_out1.img"
arcpy.env.workspace = ws2
ras_names = arcpy.ListRasters()
for ras_name in ras_names:
    name, ext = os.path.splitext(ras_name)
    ras = arcpy.Raster(os.path.join(ws2, ras_name))
    mean =arcpy.GetRasterProperties_management(ras, "MEAN")
    STD=arcpy.GetRasterProperties_management(ras, "STD")
    out_ras = ras - (mean +  3 * STD )
    out_ras.save(v_name_sc_img.format(name))


texts: Defining multiple projections inside a file gdb [error 000289 failed to alter spatial reference]?

texts:
I am importing multiple shapefiles into a filegdb in ArcGIS 10, and using the following code to define the projection of the feature once imported, to the .prj of the original .shp file.

Here is my code:

for spatial in SpatialRefList:
    del gp
    import arcpy as gp
    print 'DEFINING PROJECTION'

    print 'FEATURECLASS = ' + str(featureclass)
    print 'PROJECTION = ' + str(prjfile)
    time.sleep(1)
    try:

        gp.DefineProjection_management(featureclass, prjfile)

    except:
        traceback.print_exc()
        time.sleep(2)

However the code is throwing up error 000289 - failed to alter spatial reference.


texts: I wonder if I can add SLD style to WMS or WFS layer existing in GeoServer, using OpenLayers 3?


texts:
I ran into an issue and had to reinstall QGIS 2.18, but I can't seem to get it to properly work though. Whenever I install it it says "QGIS 2.18 is already installed on your computer..." even though I delete every possible reference to QGIS in the registry, my files, etc... I continuously get this error at startup:

    Couldn't load QGIS utils.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 21, in
    standard_library.install_aliases()
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\future\standard_library\__init__.py", line 457, in install_aliases
    __import__(oldmodname)
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\future\backports\__init__.py", line 17, in
    from .misc import (ceil,
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\future\backports\misc.py", line 24, in
    from socket import getaddrinfo, SOCK_STREAM, error, socket
  File "C:\Python27\ArcGIS10.4\Lib\socket.py", line 47, in
    import _socket
ImportError: DLL load failed: %1 is not a valid Win32 application.

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

QGIS version:
2.18.3 'Las Palmas', 77b8c3d

Python path:
['C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python', u'C:/Users/ac/.qgis2/python', u'C:/Users/ac/.qgis2/python/plugins', 'C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\nose-1.3.3-py2.7.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\certifi-14.05.14-py2.7.egg', 'C:\\Python27\\ArcGIS10.4\\Lib', 'C:\\Python27\\ArcGIS10.4\\DLLs', 'C:\\Python27\\ArcGIS10.4\\Lib\\lib-tk', 'C:\\another-library', 'C:\\PROGRA~1\\QGIS2~1.18\\bin\\python27.zip', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGIS2~1.18\\bin', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGIS2~1.18\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg']

I've tried suggestions on similar posts but they're all from older versions of QGIS and are no longer relevant. I'm not sure if this has anything to do with it but I do also have Arc installed and I use that primarily. I know Arc installs Python a bit differently since it's native.


texts: Issues with ArcGIS JSAPI in a jQuery UI Accordion?

texts:
Has anyone had success using an ArcGIS Javascript API map within a jQuery UI accordion? I have a simple example here:

http://chadwilcomb.com/maps/test/AGS_jQueryUI.htm

There is an offset between the mousedown location and the map. For instance, if you SHIFT + Drag to zoom, you will see the zoom box is severely offset from the mouse location. Also, if you double-click to zoom, the map centers on a different location than where you double-clicked. I have checked all the CSS settings in Firebug and nothing seems to really help. Anyone have any other ideas?

Additional info:

I originally experienced this problem within an ASP.NET MVC3 View, but as you can see it still occurs in a simple HTML file.

The original map had a mouseover tooltip dialogue on a feature layer and that worked as expected.

The same behavior exists when using the Dojo Accordion as well...

http://chadwilcomb.com/maps/test/AGS_DojoAccordion.htm


texts:
I have two raster layers.
Layer A has a cell size of 1020 m. There are 31,840 rows in the attribute table, each corresponding to a different 1020 m x 1020m cell.
Layer B has a cell size of 30 m. The attribute table associated with this only has two entries - value (0,1) and count of value.

I want to combine the two layers such that I have the ID of the 1020 m x 1020 m cell from layer A that spatially overlaps the values from layer B. After reading How to find the area of raster values within polygon zones? I used the "combine" tool. But the number of rows in the output are not what I would expect. Instead of the new layer having the same number of rows as layer A (31,840), there are considerably fewer (only ~22,000). The graphic below shows an example - two similar areas with roughly the same amount of missing data in layerB. An output from the Combine tool was produced for one of these areas, but not for the other! I don't understand what could have gone wrong! Any ideas?


texts: Export a Sketchup Model to ArcGIS's Multipatch?

texts:
I have a very large model made in Google Sketchup v 7.1. I would like to export this to a Multipatch featureclass, so that I can view it in ArcScene, with other 3D data.
I have tried several things, but none of them have worked. I have tried the following things:

Used the Sketchup Plugin for ArcGIS; For this I have saved the model as a sketchup 6 model, and on a machine which has Sketchup Pro 6 & ArcGIS 9.2, used the export as ESRI Multipatch. This works for small and simple models, but fails on my Large Model.
Exported the Model as KMZ, and then tried to use the QuickImport using the InterOp Extension of ArcGIS.
Exported the Sketchup 6 Model as Google Earth 4 KMZ. Then, I tried to use the QuickImport using the InterOp Extension of ArcGIS. This does import the KMZ to Multipatch, but it does not fall where I expect it to fall. Its size is very large as well (when compared to actual or expected size).

Has anyone else managed to do something like this?