texts:
I'm trying to make a supervised classification program in GEE. But it gives me the error: Layer name is not a string: ee.Classifier

// Define a region in which to generate a sample of the input.
var region = ee.Geometry.Rectangle(-71, -1.5, -67, 2.2);

// Display the sample region.
Map.setCenter(-68.50, 0.43, 7);

//select sentinel image
var img = ee.Image('COPERNICUS/S2/20170813T145731_20170813T145912_T19NEA');

//merging samples
var amostra = clareira.merge(agua).merge(mata).merge(nuvem);
//Map.addLayer(amostra);

//selecting bands
var bandas = ['B4','B3','B2'];

//creating training set
var training = img.select(bandas).sampleRegions({
  collection:amostra,
  properties:['class'],
  scale:10
});

//selecting classification algorithm
var classificacao = ee.Classifier.cart().train({
  features: training,
  classProperty:['class'],
  inputProperties:bandas});

var result = img.select(bandas).classify(classificacao);

//Map.addLayer(img, {bands:'B4,B3,B2'});
Map.addLayer(result, {min: 0, max:3, palette:['magenta', 'blue', 'green', 'white']},classificacao);


texts: Website, tool or app to produce colored map of the United States?

texts: Does anyone know of a website, Windows program or iOS app that has a clickable map of the United States and colors in the states that are clicked on? Looking to produce a map of states I have visited or lived in. Ideally it would be nice to be able to save and print this image as well. Thanks.


texts: Is there a way to render/show ACAD drawing in ArcMap correctly?

texts:
I would like to extract only certain lines from an ACAD drawing but it seems that a lot of lines are being dropped when the dwg drawng is viewed and rendered in ArcMap. (see below)  Is it possible to view an autocad drawing in ArcMap the same way it is rendered in AutoCAD or TurboCAD.

Intersection as viewed in TurboCAD:

Same file and the same intersection viewed in ArcMAP:

In ArcMap when converted to FGDB Feature Class:


texts:
I have a number of pedestrian trajectories (green) that I would like to snap to the nearest road (brown).  The roads are downloaded from OpenStreetMap and have been dissolved by name with

CREATE TABLE dissolved AS
SELECT osm_roads.name, ST_Union(osm_roads.geom) FROM osm_roads GROUP BY osm_roads.name;

Now I am left with two tables - dissolved and trajectories

My goal is first, to snap each trajectory to the nearest road and then, find the "mileage" walked along each road by

ST_LENGTH(ST_Intersection(dissolved.geom, snapped.geom))

How would I perform the first part of the query - snapping each individual trajectory to the nearest road - with ST_Snap - and then group ST_LENGTH(ST_Intersection()) by road name?


texts:
I'm playing around with Mapnik and trying to add a PNG image over a map created from an XML file. This is my XML file:

<Map background-color="white" minimum-version="0.7.2" srs="+init=epsg:3395">
    <Style name="population">
        <Rule><PolygonSymbolizer fill="#c7e9b4" /></Rule>
    </Style>

    <Layer name="countries" srs="+init=epsg:3395">
        <StyleName>population</StyleName>
        <Datasource>
            <Parameter name="type">shape</Parameter>
            <Parameter name="file">world_merc</Parameter>
        </Datasource>
    </Layer>
</Map>

And my code is like this:

import mapnik

m = mapnik.Map(600,400)
xml = 'population.xml'
mapnik.load_map(m, xml)

s = mapnik.Style()
rule = mapnik.Rule()
sym = mapnik.RasterSymbolizer()
rule.symbols.append(sym)
s.rules.append(rule)
m.append_style('Raster Style', s)

lyr = mapnik.Layer('GDAL Layer')
#lyr.srs = '+init=epsg:3395'
lyr.datasource = mapnik.Gdal(file='data.png', bbox=(0,0,1000,1000))
lyr.styles.append('Raster Style')
m.layers.append(lyr)

m.zoom_all()
mapnik.render_to_file(m,'map.png')

When that lyr.srs line is uncommented, the image is not shown, only the XML layer is there. When it is commented out, only a small part of the image is shown above the XML layer, the lower left corner of the image seems to be at the position of the origin of the SRS.

I just think that PNG files don't have info about where they should be, so the parameter bbox of the Gdal function seems to be good; however, it appears that that parameter doesn't have any effect in my code. The documentation is for mapnik 2.2 and my mapnik version is 3.0 so I'm not sure though.

Could someone tell me whether the parameter bbox has the same meaning as I think, and how I can control the position and size of the image displayed on the map?


texts:
I am working with ACS 5-Year Data at the block group level in ArcGIS, and was hoping someone could explain what one of the elements on my attribute table means.

I know that the B00001e1 is some type of population estimate, but I am confused about whether or not it is a population estimate change, or population estimate totals for each block group.


texts: Displaying real-time CSV data in OpenLayers 3?

texts:
is there a good way to display real-time csv data with OpenLayers 3?

http://geofon.gfz-potsdam.de/eqinfo/list.php?datemin=&datemax=&latmin=&latmax=&lonmin=&lonmax=&magmin=&fmt=html&nmax=10&fmt=csv

ol2 had the opportunity to display csv data, but how do i do that in ol3?


texts: Where can i download a good World image in GeoTIFF?

texts:
I can either be from MERIS (ESA) or from MODIS (Nasa). I already reached these sites:

http://www.unearthedoutdoors.net/global_data/true_marble/download

http://www.eosnap.com/public/mosri/2009/world/world-mosaic-final-100.jpg

But from the first I was not able to download, and from the second it is not available in tiff.

Thanks


texts: How to customize DB Manager default parameters

texts:
I'm using QGIS 2.10 and DB Manager plugin.

I regularly use this plugin to display tables from a PostgreSQL - PostGIS database. I particularly use the integrated SQL Window to filter a table's content with a SQL query. Then, in order to display the query's result in the QGIS map canvas, it is asked to specify Column with unique integer value and Geometry column (see screenshot below). These parameters are automatically filled with a column name and I have to change the values to gid and the_geom each time I use this...

Is there a way to specify default values for these parameters, as gid and the_geom are always the columns to be used for every table in my database?


texts: cesium-terrain-builder handle .tif to .terrain not correctly

texts:
I have a heightmap.tif on system WGS84 file in C:/Users/ctb/data/

Tilesets dir: C:/Users/ctb/data/tilesets/

Then run docker CLI: docker run -v /c/Users/ctb/data:/data -t -i homme/cesium-terrain-builder:latest bash

and run command: ctb-tile -o /data/tilesets/terrain/uo_wgs84 /data/heightmap.tif

In cesium project I get this:

but this scale is different every 5 than what should be.

Then I try in Mercator: gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 heightmap.tif heightmap2.tif and build again.

For compare heightmap.tif open in "MicroDEM":

Is it possible to do something in this situation?

EDIT:

GdalInfo for heightmap.tif:

root@9b66b4a58f72:/# gdalinfo -stats /data/heightmap.tif
Driver: GTiff/GeoTIFF
Files: /data/heightmap.tif
       /data/heightmap.tif.aux.xml
Size is 4096, 3430
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (82.888233019881397,55.458638096454301)
Pixel Size = (0.000011216234429,-0.000006375170268)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  82.8882330,  55.4586381) ( 82d53'17.64"E, 55d27'31.10"N)
Lower Left  (  82.8882330,  55.4367713) ( 82d53'17.64"E, 55d26'12.38"N)
Upper Right (  82.9341747,  55.4586381) ( 82d56' 3.03"E, 55d27'31.10"N)
Lower Right (  82.9341747,  55.4367713) ( 82d56' 3.03"E, 55d26'12.38"N)
Center      (  82.9112039,  55.4477047) ( 82d54'40.33"E, 55d26'51.74"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Min=56.131 Max=98.161
  Minimum=56.131, Maximum=98.161, Mean=63.756, StdDev=6.603
  NoData Value=-32767
  Overviews: 2048x1715, 1024x858, 512x429, 256x215
  Metadata:
    STATISTICS_MAXIMUM=98.160507202148
    STATISTICS_MEAN=63.756109311469
    STATISTICS_MINIMUM=56.131080627441
    STATISTICS_STDDEV=6.6032626523927

P.S. I am sorry for my knowledge of the English language.


texts: GeoServer faling to connect to web admin page or homepage

texts:
I'm trying to get started with learning GeoServer but can't get it to work.  I downloaded 2.6.1 for windows 7.  Install runs fine and my Java is updated.  When I tried to connect to either the web admin page or the homepage I get an "Application not Found" error.  I've tried setting it to different ports but same thing.

Shot from GeoServer Starter command window


texts: Missing something when using custom CRS to define projection for NOAA QPF

texts:
I'm trying to display NOAA QPF data and this: https://spatialreference.org/ref/sr-org/6770/
suggests the proj string should be:

Proj4: +proj=longlat +a=6371200 +b=6371200 +no_defs

which matches:
https://www.wpc.ncep.noaa.gov/html/mapprojection.txt

GEOGCS["GCS_Sphere_EMEP",
    DATUM["D_Sphere_EMEP",
        SPHEROID["Sphere_EMEP",6371200.0,0.0]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]

However, I get this as output:

Obviously, I'm doing something wrong. Can anyone point me in the right direction?  My overall goal is to ingest the WFS from here: https://idpgis.ncep.noaa.gov/arcgis/services/NWS_Forecasts_Guidance_Warnings/wpc_qpf/MapServer/WFSServer  into QGIS. Thanks in advance!


texts:
How does one subclass QgsMapToolDigitizeFeature?
I'm trying to create a map tool in QGIS using Python. The tool will need to draw a polygon (and points and lines) in the same manner that QGIS does natively (including obeying any topological rules). Once the feature is placed I will then provide any attributes necessary for the feature.
I've created a tool for adding points using QgsMapToolEmitPoint but can't help but feel that I'm doing a lot of work that is repeated within QGIS its self.
So what is the minimum code required to subclass QgsMapToolDigitizeFeature that allows a user to draw a polygon. Just doing:
class TestPolygonTool(QgsMapToolDigitizeFeature):

    def __init__(self, canvas, cdw):
        self.canvas = canvas
        self.cdw = cdw
        super(TestPolygonTool, self).__init__(self.canvas, self.cdw, QgsMapToolCapture.CaptureNone)

    def activate(self):
        super(TestPolygonTool, self).activate()

and then activating via:
testTool = TestPolygonTool(iface.mapCanvas(), iface.cadDockWidget())
iface.mapCanvas().setMapTool(testTool)

does not appear to be enough.
And more generally where did you find this information? Nobody seems to talk about doing this.


texts: I want to cut polygons in a polygon layer.

texts:
I want to programmatically cut polygons already loaded in arcmap. For instance if originally I have 1 polygon feature and 1 row in attribute table. When I cut them, I should see 2 independent features and 2 rows in attribute table. I tried using the snippet from the reference with no luck

Any suggestions how i can get this using C# ArcObjects.
Any suggestions as to how i can achieve this. It


texts:
My script is terribly underdeveloped, I am trying to select cities within a country and then order these cities by their "POP_RANK" = population rank. I am not sure whether to use a searchcursor or just select the features and then use a sort function.

import arcpy
from arcpy import mapping
#set workspace
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
#set the dataframe
fc = ("N:\Lab13\cities.shp")
cities = arcpy.mapping.ListLayers(mxd, "Cities")[0]
cities.showLabels = True
def countryname():
    country = arcpy.GetParameterAsText(0)
    country_layer = arcpy.MakeFeatureLayer_management(fc, "country_lyr")
    arcpy.SelectLayerByAttribute_management(cities, "NEW_SELECTION", "CNTRY_NAME = '{}' ".format(country))
    arcpy.mapping.ListDataFrames(mxd)[0].zoomToSelectedFeatures()
    arcpy.RefreshActiveView()
    arcpy.Delete_management("country_lyr")
try:
    in_dataset = "CNTRY_NAME"
    out_dataset = "country_sort"
    sort_fields = [["POP_RANK", "ASCENDING"]]
    arcpy.Sort_management(in_dataset, out_dataset, sort_fields)
    arcpy.AddMessage("IT WORKED")


texts: Why does osr.SpatialReference return different values in python2 and 3?


texts: Land Cover Classification Sentinel-2 deep learning model

texts:
I want to create a land cover using the model in the link (https://www.arcgis.com/home/item.html?id=afd124844ba84da69c2c533d4af10a58). For this you need to mosaic Sentinel-2 images using the tool in the link (https://www.arcgis.com/home/item.html?id=e6e1f20cb0374d28a6eed24f5c2ff51b). When I run the land cover model using the mosaic, almost all land is classified as inland waters. I get more accurate results when I do the images without mosaic, just using XML file in Sentinel-2 image folder. However, the land cover model uses 60 m resolution images. The deep learning model is designed to use 10m resolution multispectral images by mosaic.
Has anyone used the model and experienced this problem before?


texts: Removing extra spaces from Attribute Assistant calculation

texts:
Using the Attribute Assistant in ArcGIS 10.5, I auto-populate a field called FULLADDR by entering values for the fields below:

([ADDRNUM]+" "+[ROADPREDIR]+" "+[FULLNAME]+" "+[ROADPOSTDIR]+" "+[UNITTYPE]+" "+[UNITID])

If there is no ROADPREDIR, a space is entered into the FULLADDR field where the ROADPREDIR would be if it existed.  I would like to ensure that if there is no ROADPREDIR value there will be no space between fields ADDRNUM and FULLNAME from the equation above.


texts: os module can't find time files were created, because "system can't find file specified"

texts:
I have ~200 layers in a Map Document, and I am trying to find the dates that these layers were created on (so I can export to database later). Some of these layers are housed in groups and sub-groups and also in gdbs, and I think this is creating a problem when I try to extract the metadata.

# import the modules.
import arcpy
import os
import time

# specify path of map document where layers are located.
path = 'Location/1_STATEWIDE/'
er = path + "Tool 10_2.mxd"

# list layers using mapping module. I used this for earlier versions of the code where I had to extract layer longnames and datasources.

mxd1 = arcpy.mapping.MapDocument(er)
layers1 = arcpy.mapping.ListLayers(mxd1)

# loop asks that for every layer in the Map Document, if the layer has a datasource, the os module will proceed to get creation time of the layer.
The datasource condition exists so that group layers and sub-group layers will not give any results or errors; I just need times for the layers.

for lyr in layers1:
    if lyr.supports("DATASOURCE"):
        sources = lyr.dataSource
        # need path to call time, so print list of paths before time commands.
        print os.path.getctime(sources)

This code gives the following error:

WindowsError: [Error 2] The system cannot find the file specified:     u'\\\\Location\\Terrapins_horseshoecrabs\\2002terrapindata\\terrapin_survey2002.gdb\\beach_survey_2002'

So this must be because I am trying to find the creation time for the files in a geodatabase. How can I correct this code to avoid this error?


texts: Copy column from Excel to attribute table

texts: Is there any way to copy a column from excel to an attribute table's column in ArcGIS?


texts: Cause of ExecuteError: ERROR 010067 from arcpy.sa.CostPath?

texts:
By using the code below that creates least-cost paths between several polygons, I obtain this error message:

ExecuteError: ERROR 010067: Error in executing grid expression. Initalization failed

when I execute this code line:

outCostPath = CostPath("Selected_patch_end", "outCostDist.tif", "outCostLink.tif", "EACH_ZONE","FID")

Sometimes, when I re-execute my script, it works.

However, when the code works, I always obtain strange files and folders every time that I execute the same code line above:

Here is an insight of my dataset:

Here is my code:

 import arcpy
 from arcpy import env
 from arcpy.sa import *
 arcpy.CheckOutExtension("Spatial")
 arcpy.env.overwriteOutput = True
 arcpy.env.workspace = r"C:\Users\LCP"
 arcpy.env.extent = "costsn.tif"

with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInPatches_start:

   for rowStart in rowsInPatches_start:

      ID_patch_start = int(rowStart[3])

      ## Define SQL expression for the fonction Select Layer By Attribute
      expressionForSelectInPatches_start = "GRIDCODE=%s" % (ID_patch_start)

      ## Process: Select Layer By Attribute
      arcpy.MakeFeatureLayer_management("selected_polygons.shp", "Selected_patch_start", expressionForSelectInPatches_start)

      ## Process: Cost Distance
      outCostDist=CostDistance("Selected_patch_start", "costsn.tif", "", "outCostLink.tif")

      ## Save the output
      outCostDist.save("outCostDist.tif")

      with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInSelectedPatches_end:

         for rowEnd in rowsInSelectedPatches_end:

            ID_patch_end = int(rowEnd[3])

            ## Define SQL expression for the fonction Select Layer By Attribute
            expressionForSelectInPatches_end = "GRIDCODE=%s" % (ID_patch_end)

            ## Process: Select Layer By Attribute in Patches_end
            arcpy.MakeFeatureLayer_management("selected_polygons.shp", "Selected_patch_end", expressionForSelectInPatches_end)

            ## Process: Cost Path
            outCostPath = CostPath("Selected_patch_end", "outCostDist.tif", "outCostLink.tif", "EACH_ZONE","FID")

            ## Save the output
            outCostPath.save('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".tif")

            ## Write in a file .txt
            outfile=open('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".txt", "w")
            rowsTxt = arcpy.SearchCursor('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".tif")
            for rowTxt in rowsTxt:
               value = rowTxt.getValue("Value")
               count = rowTxt.getValue("Count")
               pathcost = rowTxt.getValue("PATHCOST")
               startrow = rowTxt.getValue("STARTROW")
               startcol = rowTxt.getValue("STARTCOL")
               print value, count, pathcost, startrow, startcol
               outfile.write(str(value) + " " + str(count) + " " + str(pathcost) + " " + str(startrow) + " " + str(startcol) + "\n")
            outfile.close()


texts: Reference materials for learning Google API V3 with PHP and MySQL (tracking)

texts:
I am new to this Google API. Now i wanted to learn Google API with PHP and MySQL.

Actually i wanted to create an application  that will track x,y coordinates coordinates on real time bases for android mobile. There will be many Android mobiles which will send data on real time bases and using php and mySQL i wanted to display it on Google map and refresh display and data after certain time interval.

So with respect to this tracking system part can any one suggest me any references material and web site.


texts: Clipping interpolation layer to a coastal boundary in QGIS 3.6?

texts: I have created an interpolation layer and want to clip it so I don't see the area outside the coastal boundary. map showing area to be clipped:


texts: Mass-update the projection in compressed GeoTIFF files

texts:
I have a huge set of compressed GeoTIFF files. Someone forgot to set the projection while creating the GeoTIFF files and I would like to add it now.

If I use

  gdalwarp -t_srs source.tif target.tif

the result is an uncompressed file. I would like to avoid the (de)compressing if possible, since I expect it will degrade the final result due to the recompression and because it will require more time to decompress and compress each file again.

I assume that the projection information is stored in the header.

Is there an easy way to just update the projection information in a GeoTIFF file?


texts: Finding points of sudden large changes in stream gradient QGIS

texts:
In QGIS I am trying to divide a streams shapefile at points at where there are sudden large increases or decreases in gradient.

For example at the red circle the gradient of the stream increases from about 3% to nearly 15%. I would like to split the stream feature at sudden changes like these.

I have gone through and done this by hand with many streams using the editor split features tool but this is slow. I am looking for a more efficient way to find these point and split the feature.


texts: Creating fields in attributes table using other layers using QGIS

texts:
I want to create a grid over my map and for each box in that grid I want to record attributes from other layers in my project; or I want to add columns in the attributes table for the grid layer where these columns contain data about other layers in my project.

For example if my grid layer is named grid and I have another layer called location which consists of points, then I want to add a column to the grid attributes table that gives me the total number of points in location that coincide with each box of grid. This is a simple example, but I need to do this for many layers to create many fields in the grid attributes table. My hope is export this attributes table for economic analysis.


texts: Can you access ArcScan from a python script?

texts:
I'm using a python script with ArcPy in ArcGIS 10, and want to access some ArcScan functions.

Is this possible? I can't see anything in the documentation about it, but I can't really find any very good ArcPy documentation.

Update: Just to add some extra information: I want to run the centreline vectorisation method (see here) from an ArcPy script. Is that possible?


texts:
For the first time, I am using NetCDF data files containing daily data about soil moisture at a global level.

I want to be able to obtain monthly mean  value of soil moisture for each grid cell.

There is one NetCDF file for each day, and each file contains four bands. I only need the first band which contains the data I want.

Using QGIS 1.8, I encounter the following problems:

1. Loading the rasters

When I want to load the 31 rasters related to January (example), I can select the 31 files, but then, I have to select the band I am interested in for each layer, one by one. Is there a way to load all of them, but automatically select the band I want?

2. Converting the layers

The data comes in a floating number format. I detected some issues when it comes to sum values. However, when I convert the layer to a .tiff raster, it seems to transform the value into an integer. This is what I want, but GDAL doesn't want to batch convert all the files into .tiff, either because of the file names of because of the multiband issue. Is there a way to batch convert NetCDF files into .tiff?

3. Batch summing the layers

Finally, even if I have the layer in .tiff format, I wonder if there is an easier way to sum all the values than by using RasterCalculator. Is there a tool than batch sums up rasters within a folder?
Also, I noticed that I cannot use RasterCalculator directly on the NetCDF layers because of the file names. If the solution is on this trail, is there a way to batch change the name of a layer within a NetCDF file?


texts: How to resolve GRASS error "Execution failed: 'gproj -p' "?

texts: I installed the WinGRASS-6.4.0-1-Setup.exe 55M from http://grass.osgeo.org/grass64/binary/mswindows/native/. I'm using a x64. After installation, I tried opening the demoLocation but it failed and displayed the error: Execution failed: 'gproj -p'. How do I resolve it?


texts:
I'm working on a geoserver instance. A style directory can be found in two places:

A) /styles
B) /data/workspaces/default/styles

Both contain the same set of styles except the one in the root contains the generic line, point, etc. styles.
My guess is the latter exists for portability reasons? My questions are:

which one is used? Or how does Geoserver decide which style is picked if styles are identical in both directories?
To tidy things up, in which folder should I remove styles to only have one copy. For example Keep generic in A) and the others in B) or only use A) and delete the style folder in workspaces  (B) ?


texts: AutoCAD Map 3D does not show panels and icons?

texts: I have some problems with autocad 2014, 3D map, this not show me de icons on the top, I use windows 10.


texts: Terrain correction / Orthorectification with gdal

texts:
Since MODIS Imagery Orthorectification (Open Source) had no answer, I will try to get one here with a similar problem:

I got some raster datasets with GCPs in lat/lon, but without height information. When I use gdal to project/re-project the data, I get an offset of approx. 2km, which I think comes from terrain distortion.
Using ESA's SNAP Toolbox, I get the results I want, but they are using the actual SRTM DEM in the background.

Is there a way to provide gdal with the raster, the GCPs and a DEM to get everything done "at once" (or in several steps)?

So far, my approach (which gives the offset) looks like this:

gdal_translate -gcp pix1 line1 x1 y1 -gcp pix2 line2 x2 y2 (...) input.tif output.tif

The data are SAR raster from ERS-1, which might have something to do with the offset due to the viewing angle of the sensor


texts: QGIS NetCDF raster and shape vector not overlapping in EPSG: 27700

texts:
I have got a raster file and a vector file of a UK which I try to overlap to use zonal statistics. The files according to the data providers are projected in British National Grid (EPSG 27700). However, when I load the 2 files, QGIS first fails to recongize the correct CRS and creates a custom CRS for each.

For the vector:
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601272 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

For the raster:
+proj=tmerc +lat_0=49 +lon_0=0 +k=1 +x_0=400000 +y_0=-100000 +ellps=airy +units=km +no_defs

The resulting image I get

Although the maps are overlapping visually, zonal statistics does not output anything. When I set both files CRS to EPSG 27700, I get the following:

The raster file is the small point in the centre, it is much smaller and is rotated.

How could I align the two maps?


texts: Python syntax returning only part of the if/elif statement in Field Calculator?

texts:
Using ArcGIS 10.7 desktop.
I want to reclassify a land use field using Field Calculator. Grasslands are conditionally reclassified using a 'veg_type' field (Avena, or Bromus) and a CWHR land use field (grassland values = Annual Grassland, AGS, Perennial Grassland). All other values in the CWHR field should stay the same.
When I run the Python Parser syntax the Grassland_Reclass field returns "Avena Grassland" (gridcode = 1) for the assigned land use values, but also for grasslands categorized as "Bromus' (gridcode = 0). The default CWHR values are added, but the elif condition (gridcode = 0) is not returned. I've also tried with the veg_type field but I get the same result.
def Reclass(gridcode, CWHR_TYPE):
  if gridcode == 1 and CWHR_TYPE == "Annual Grassland" or CWHR_TYPE == "AGS" or CWHR_TYPE == "Perennial Grassland":
    return "Avena Grassland"
  elif gridcode == 0 and CWHR_TYPE == "Annual Grassland" or CWHR_TYPE == "AGS" or CWHR_TYPE == "Perennial Grassland":
    return "Bromus Grassland"
  else:
    return CWHR_TYPE

Expression: Grassland_Reclass =
Reclass(!gridcode!, !CWHR_TYPE!)


texts:
I am currently trying to generate an empty database in pgAdmin 4.29 using ili2pg. The Interlis model SIA405_ABWASSER_2015_LV95 is available here. The ili2pg-command I used is
java -jar ili2pg-4.4.4.jar --schemaimport --dbhost myhost --dbusr postgres --dbpwd ******
--dbdatabase mydb --dbschema test --setupPgExt --coalesceCatalogueRef --createEnumTabs
--sqlEnableNull --createFkIdx --coalesceMultiSurface --coalesceMultiLine --coalesceMultiPoint
--coalesceArray --beautifyEnumDispName --createGeomIdx --createMetaInfo --expandMultilingual
--createTypeConstraint --createEnumTabsWithId --createTidCol --importTid --smart2Inheritance
--createBasketCol --defaultSrsCode 2056 --trace --models SIA405_ABWASSER_2015_LV95 -modeldir my_dir

However, this command generates a java.lang.NullPointerException. According to the trace log, the error stems from an update in T_ILI2DB_SETTINGS after the EnumValueMap is done:
Info: createEnumValueMap(): SELECT iliCode,T_Id FROM test.sia4055_lv95sia405_bwssr_lk_haltung_lagebestimmung WHERE thisClass = 'SIA405_ABWASSER_2015_LV95.SIA405_Abwasser_LK.Haltung.Lagebestimmung' (EnumValueMap.java:52)
Info: updateMultiEnumTableWithId(): INSERT INTO test.sia4055_lv95sia405_bwssr_lk_haltung_lagebestimmung (seq,iliCode,itfCode,dispName,inactive,description,thisClass,baseClass) VALUES (?,?,?,?,?,?,?,?) (TransferFromIli.java:1504)
Info: readColumnTab(): SELECT tablename,subtype,columnname,tag,setting FROM test.T_ILI2DB_COLUMN_PROP (DbExtMetaInfo.java:144)
Info: saveColumnTab(): INSERT INTO test.T_ILI2DB_COLUMN_PROP (tablename,subtype,columnname,tag,setting) VALUES (?,?,?,?,?) (DbExtMetaInfo.java:91)
Info: readTableTab(): SELECT tablename,tag,setting FROM test.T_ILI2DB_TABLE_PROP (DbExtMetaInfo.java:232)
Info: saveTableTab(): INSERT INTO test.T_ILI2DB_TABLE_PROP (tablename,tag,setting) VALUES (?,?,?) (DbExtMetaInfo.java:184)
Info: readIliFiles(): SELECT filename,iliversion,modelName FROM test.T_ILI2DB_MODEL (TransferFromIli.java:457)
Info: addModels(): INSERT INTO test.T_ILI2DB_MODEL (filename,iliversion,modelName,content,importDate) VALUES (?,?,?,?,?) (TransferFromIli.java:585)
Info: updateSettings(): INSERT INTO test.T_ILI2DB_SETTINGS (tag,setting) VALUES (?,?) (TransferFromIli.java:736)
java.lang.NullPointerException
    ch.ehi.ili2db.metaattr.MetaAttrUtility.visitElement(MetaAttrUtility.java:187)
    ch.ehi.ili2db.metaattr.MetaAttrUtility.visitElement(MetaAttrUtility.java:194)
    ch.ehi.ili2db.metaattr.MetaAttrUtility.visitElement(MetaAttrUtility.java:194)
    ch.ehi.ili2db.metaattr.MetaAttrUtility.visitElement(MetaAttrUtility.java:194)
    ch.ehi.ili2db.metaattr.MetaAttrUtility.updateMetaAttributesTable(MetaAttrUtility.java:160)
    ch.ehi.ili2db.base.Ili2db.runSchemaImport(Ili2db.java:1291)
    ch.ehi.ili2db.base.Ili2db.run(Ili2db.java:227)
    ch.ehi.ili2db.AbstractMain.domain(AbstractMain.java:554)
    ch.ehi.ili2pg.PgMain.main(PgMain.java:70)
Finished with errors!

using ili2gpkg results in the same error.


texts: Preventing reflective surface in rasterVis

texts:
I have been experimenting with digital terrain model files in ESRI ASCII grid raster format. Despite zero experience with this type of data I found it very easy to load into R, as per code below. The plot3D function from the rasterVis package defaults to a nice looking plot but the model is rather reflective. Given that it's supposed to be a landscape this shiny surface is not really appropriate.

I am looking for a more natural (if such a word can be used about a 3D rendering of a set of figures), matt, topological plot.

I imagine there must be ways around this but I have no experience using raster and seek some pointers.

How can I reduce or remove the reflectivity of the plot?

library(raster)
library(rgdal)
library(rasterVis)

foo <- raster(readGDAL("my.dtm.asc"))
plot3D(foo)


texts: Querying attachment info from ArcGIS map service?

texts:
Like feature layers have the queryAttachmentInfos() method to fetch attachment info, is there something like this to fetch attachment info from a map service (ArcGISDynamicMapServiceLayer)?

My app flow is as follows

query the map service based on some user defined criteria
set info templates for the returned features (attributes + geometry) (this is also where i specify showAttachments: true but again no luck)
add features to map.graphics and show the info window (custom popup in my case)

I have also tried using the attachment editor and the attribute inspector but with no luck.

Should i use the relationship query task? The problem with using the relationship query, i assume, would be that it would be impossible to know if the related record is an attachment or some other related record.


texts: How to add other images as TIFF overviews?

texts:
TIFF files can contain overviews for quicker loading/displaying at small scales (=more "zoomed out").

These can be generated and added with GDAL's gdaladdo command. gdaladdo will compute the overviews by summarizing blocks of pixels as individual pixels. (Thus why the overviews are also sometimes called "pyramids".) This makes sense for orthophotography or for scalar-raster-geodata-encoded-as-image TIFFs (e.g. a DEM), but if the TIFF is a map, I think it might be more suitable to use smaller-scale (=more "zoomed out") maps (if available) as overviews that has e.g. a lower level of detail, omits some features and has captions that are still readable at that display size.

Is there a way to use separate images to create the overviews to be included in a TIFF, maybe even with GDAL command line tools?

Would doing so break any overview use-cases or tools that are also relevant for maps?


texts:
I'm trying to build a model using Model Builder in ArcGIS 10.1 that starts off by joining two shapefiles together based on a common field. I would like to make this open ended (not hardcoded) so I can use it repeatedly.

Problem is when I go to select my shapefiles when I run it, I can add one, but not both (I only get an option to select "data elements").

Likewise, when I try to edit the model and insert or hardcode the shapefiles into it, I get an error message at Add Join that says

ERROR 000840
The value is not a Raster Layer. ERROR 000840: The value is not a Raster Catalog Layer. ERROR 000840: The value is not a Mosaic Layer.

When I go through this same workflow in ArcMap I get the results I am looking for.

Any idea what is causing these errors, how to correct them or a work around?

So we are clear, the following workflow is what I'm trying to automate:

1.  Join [CURRENTMONTH] shapefile to [PREVIOUSMONTH] shapefile.
2.  Perform a select by attributes:
    a.  Select From “[PREVIOUSMONTH]_[CURRENTMONTH]“ WHERE “[CURRENTMONTH].SEGMENT_ID” IS NULL
3.  Remove Join.
4.  Export Selected Features to new shapefile (deletes.shp)

The second half of the model reverses the process and attempts to accomplish the following:

1.  Join [PREVIOUSMONTH] shapefile to [CURRENTMONTH] shapefile.
2.  Perform a select by attributes:
    a.  Select From “[CURRENTMONTH]_[PREVIOUSMONTH]” WHERE  “[PREVIOUSMONTH].SEGMENT_ID” IS NULL
3.  Remove Join.
4.  Export Selected Features to new shapefile (added.shp)

UPDATED TO SHOW MODEL DIALOG

Updated to show new model


texts:
As asked in previous question: Displaying icon name next to icon using OpenLayers? , now I need to change data source to the JSON API data, and used jquery plugin bypass CORS (http://www.ajax-cross-origin.com/install.html).
The code part i changed from previous question jsfiddle:

var apiURL = 'https://map.navifleet.pl/api.php?api_key=[api-key]&type=1';

$.ajaxSetup({
    crossOrigin: true,
    proxy: 'proxy.php'
});

$.getJSON(apiURL, null, function(data){
        console.log(data);
    console.log(JSON.stringify(data));
  var transform = ol.proj.getTransform('EPSG:4326', 'EPSG:3857');
  data.forEach(function(item) {
     console.log(item);
     console.log(JSON.stringify(item));
    var feature = new ol.Feature(item);
    feature.set('name', item.name);
    var coordinate = transform([parseFloat(item.lon), parseFloat(item.lat)]);
    var geometry = new ol.geom.Point(coordinate);
    feature.setGeometry(geometry);
    mapSource.addFeature(feature);
  });
});

The code returns Uncaught TypeError: data.forEach is not a function.
Updated code: code

Example JSON Data from API:

{
    "success": true,
    "data": [{
        "device": {
            "id": 1555,
            "name": "Iveco Daily 2015",
            "icon_color": "009900",
            "mark": "IVECO",
            "model": "DAILY",
            "license": "K1 NAVI",
            "mileage": 20966.707
        },
        "data": [{
            "gps_timestamp": 1542282655,
            "lon": 18.5842,
            "lat": 49.7232,
            "speed": 101,
            "direction": 239,
            "distance": 0.8741,
            "total_distance": 199918
        }]
    }, {
        "device": {
            "id": 7102,
            "name": "Iveco Daily 2018",
            "icon_color": "1873BA",
            "mark": "Iveco",
            "model": "Daily",
            "license": "",
            "mileage": 0
        },
        "data": [{
            "gps_timestamp": 1541508344,
            "lon": 20.0676,
            "lat": 50.0865,
            "speed": 0,
            "direction": 258,
            "distance": 0,
            "total_distance": 322.622
        }]
    }]
}


texts: How to Reduce Decimal Precision of Image Value in ArcGIS Desktop?

texts:
I'd like to ask about the value on following image.

This is an interpolated precipitation layer using IDW method in ArcGIS Desktop. On this layer, the values is like this:

This is a raster image. After exporting the above layer into raster image, the value of the raster image has quite long decimals, specially the minimum value looks like this 4.43193612012e007.

I have not only this image but also several images with similar value. When I calculate any equation on those images, output value from that calculation cannot be accepted, they were totally wrong.

So how can I reduce the decimal precision before or after exporting to raster?


texts: Converting DWG to shapefile - include annotations in the point attribute table

texts:
I am working with a GPS system (Magnetfield system) that cannot easily export the data collected into a shapefile - it creates an individual shapefile for each point taken rather than exporting the entire job.
So, we usually export to DWG format. This creates a headache on the GIS side of things. When I try to convert the DWG to shapefile - the points collected are in one layer and the names of the points and point numbers are in the annotation layer.
Is there a way to easily add the codes and point numbers into the point layer attribute table. I don't want to have to go into the data - identify each point and create a new field and fill in each and every one, this is a big job with lots of points collected.
This seemed a lot easier when I was working with Trimble units.
I have messed around with our GPS unit and there doesn't seem to be a simple way to just export the entire job to shapefile format.


texts: QGIS 3.6 and QGIS 3.10 report different geometry types for the same shapefile

texts:
I used QGIS 3.6.3 to create a polyline shapefile, and the program reported the polyline as LineStringZM. When I opened the same shape file in QGIS 3.10.3-2, the polyline was reported as MultiLineStringZM. (The geometry type can be found by right-clicking the shapefile on the Layers panel, then "Properties", then "Information", then under "Geometry")

I am confused why the same polyline can be LinestringZM or MultiLinestringZM when opened in different versions of program?  Did I do something wrong when creating the polyline?  This is important to me, because I use python to process the polyline and different geometry type would require different coding. My code is able to process the shapefile in QGIS 3.6.3, but gags when used in QGIS 3.10.


texts:
I'm trying to rotate a box around its southwest corner using PostGIS.  However, I am not getting the expected results when doing the rotation.  I'm using SRID 3785.

The image below shows what distortion I'm talking about.

Basically, the original is a square, but the same polygon rotated 90 degrees is distorted.  Rotated by 180 is back to normal, and rotate by 270 is distorted the same way.

They should all be squares. What's going on?


texts: Setting label value for varied symbology by attribute using size as variation

texts:
In using ArcGIS Pro 2.4.1, I wanted to have symbology that represented two attributes at one time. This is possible using a variety of symbology setups, specifically, I did the following:
Created a Primary Symbology using Unique Values as shown below:

Then I changed the Vary symbology by attribute tab and used a custom expression to set the unique values to be different sizes based on my expression as shown below:

My symbology of the custom expression has resulted in the below:

How can I customize the label values from being 6, 14, 22, & 30?


texts: Join multiple rows to a single feature in QGIS

texts:
I have this kind of csv attribute table:

          ID INVENT_KPV KASVUKOHA_KOOD RINDE_KOOD PUULIIGI_KOOD VANUS
297945203711 2009.08.31             MO          1            KS    60
297945203711 2009.08.31             MO          2            KU    50
297945303711 2009.08.31             SS          1            MA    60
297945303711 2009.08.31             SS          1            KS    60
297945403711 2009.08.31             JO          1            HB    NA

As you can see, there are 3 different ID codes. I have also a map in QGIS and I want to join this attribute table and the map together but the problem is that when I use the button "info" in QGIS and click on the map then it shows me only information about one row not all of them (which have the same IDs).

For example, if I click somewhere in the map where's a patch with ID 297945303711 then it only shows me: 2009.08.31, SS, 1, MA, 60 and not 2009.08.31, SS, 1, KS, 60.


texts: How to find the image ID participating in the mosaicing in Google Earth Engine

texts:
I want to find the image ID participating in the mosaicing, here is my mosaic process. Just the partial images appearing in the mosaicing result, not all those meet the filter conditions. The color difference of the whole image is too big, I want to match each image with histogram and then mosaic it, or maybe you have a better way?
var myimage = ee.ImageCollection('COPERNICUS/S2_SR')
          .filterDate('2020-05-01', '2020-08-01').sort('CLOUDY_PIXEL_PERCENTAGE',
false).mosaic();
var myimage1 = myimage.clip(dongbei);
var Viscolor= {
bands              :['B8','B4','B3'],
min                : 0.0,
max                : 10000.0,
};

Map.addLayer(myimage1, Viscolor, ' Color (843)');

"dongbei" is the shapefile I uploaded.


texts: Converting WGS_1984_Web_Mercator_Auxiliary_Sphere to a lat long raster

texts:
I have a dataset I downloaded with the following attributes shown with gdalinfo. How can I convert it to a lat long based raster file?

Driver: GTiff/GeoTIFF
Files: d2wesri.tif
Size is 1300, 900
ERROR 1: PROJ: proj_create_from_database: Open of /home/default_user/.conda/envs/user/share/proj failed
proj_create_from_database: Open of /home/default_user/.conda/envs/user/share/proj failed
Coordinate System is:
ENGCRS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
    EDATUM[""],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037507.067200001329184,19302890.717915385961533)
Pixel Size = (30826.923076923078042,-30826.923076923078042)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-20037507.067,19302890.718)
Lower Left  (-20037507.067,-8441340.051)
Upper Right (20037492.933,19302890.718)
Lower Right (20037492.933,-8441340.051)
Center      (      -7.067, 5430775.333)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  Metadata:
    SourceBandIndex=0

I would like a lat long based raster file that would look like this:

Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-179.999988908250998,83.329027777777995)
Pixel Size = (0.005208333332500,-0.005208372298073)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-179.9999889,  83.3290278) (179d59'59.96"W, 83d19'44.50"N)
Lower Left  (-179.9999889, -55.9063889) (179d59'59.96"W, 55d54'23.00"S)
Upper Right ( 180.0000110,  83.3290278) (180d 0' 0.04"E, 83d19'44.50"N)
Lower Right ( 180.0000110, -55.9063889) (180d 0' 0.04"E, 55d54'23.00"S)
Center      (   0.0000111,  13.7113195) (  0d 0' 0.04"E, 13d42'40.75"N)

I tried using gdalwarp -of GTiff -s_srs epsg:3857 -t_srs epsg:4326 d2wesri.tif d2wesri_warped.tif
to convert the file, but it fails with

ERROR 1: PROJ: proj_create_from_database: Open of /home/default_user/.conda/envs/user/share/proj failed
ERROR 1: Translating source or target SRS failed:
epsg:3857

I uploaded a sample file here as well.


texts: Virtual Layer from SQL query is extremely slow

texts:
I have two layers, one is a grid of circles and the other is all Census Block Groups (CBG) in the state of Florida. I created a Virtual Layer out of this code
SELECT CirclesFlorida.name, group_concat(cbgFl.FIPS)
FROM CirclesFlorida
LEFT JOIN cbgFl ON ST_Intersects(circlesFlorida.geometry, cbgFl.geometry)
GROUP BY circlesflorida.name

The purpose is for each circle to have a data column which lists every CBG that falls in that circle.
Is there a faster way to do this? Maybe DB Manager? Currently the query takes 20 minutes to run and another 5 minutes to just open the attribute table of the Virtual Layer.


texts: want to burn my qgis project on cd_problem with opening pictures and changing directory

texts:
I have created a QGIS project and want to burn it on CD. For a better demonstration I have linked pictures to every point in my point-shapefiles. The problem that I am facing is, that when I open the project on a different computer via CD, everything is shown but not the pictures when I use the identify tool, as the linkage to my c-system is no longer available. How can I change the directory of my photo linkage so that it uses the folder with the pictures that is provided on the CD but no longer on the c-system.
My whole project is in one main folder that already contains the picture folder. I will burn the whole folder on the CD. Only the linkage with the pictures does not work. The rest opens up just fine. Is there any solution for that problem?


texts: Defining Parameters using Spatially Constrained Multivariate Clustering of ArcGIS Pro

texts:
I have a point shapefile with 51673 features. I am trying to cluster them using the Spatially Constrained Multivariate Clustering tool in ArcGIS Pro. After reading through the documentation, I'm confused about how to define the parameters properly to get just 5 clusters that contain around the same number of points.

I thought to put 10,000 as the minimum number per cluster, and then specify the number of constraints as 5. However I get an error when I try to run this:

  ERROR 110122: Constraint Error: Unable to solve the tool based on the
  Minimum/Maximum Constraints.  Failed to execute (Spatially Constrained
  Multivariate Clustering).

I also tried defining the maximum number as 10,000 - but this gave me 9 clusters in which the distribution of features was unequal within the groups (and none of clusters actually got close to the maximum count of 10k). Again, ideally I want 5 clusters containing a similar number of features.

How could I modify the parameters so the tool gives me this back?


texts:
iface.addVectorLayer("vector_path","vector_name","ogr")

This code does not rename the layer as "vector_name". This should not be happening, right? Is there a work around?


texts:
in an attempt to import editable data to mapbox I fail to get the data imported in a form that the system accepts.

It rejects the data with an error message stating:

  cities.geojson: Input failed. old-style crs member is not recommended
  on line 1.

When reading that I assumed somthing must have gone wrong with the conversion of my file, so I compared it to the sample file provided to MapBox without making out a difference.

Then comes the most confusing part: When trying to import the sample file I get a very similar error message from the system:

  stations.geojson: Input failed. old-style crs member is not recommended, >this object is equivalent to the default and should be removed on line 1.

Now my question is: how can import a geoJSON file into MapBox?


texts: ArcMap does not "remember" linked custom .style files

texts: I have a two sets of custom .style files that I would like to have ArcMap load every time I open the program. (Ex: the thematic color ramps from Color Brewer.) It's a small thing but I don't want to have to manually add them every time. Is there a workaround for this?


texts: Many mapfiles with few layers, or all layers in one mapfile?

texts:
I am using MapServer 6 and defining maps with mapfiles.

I have a page where 20 map compositions (~20*5=~100 layers) are possible by varing LAYER/DATA  field and some classes of layers... So I can express each map  by a mapfile (20 mapfiles) or I can put all layers in the same mapfile... This is the question,

What is better (have better performance), many mapfiles with few layers each, or one  mapfile with all (many) layers?

If the answer need some context... I am accessing Mapserver by this URL, where layer names (L1a,L2b) and mapfile name (MFc.map) can vary,

example/cgi-bin/mapserv?LAYERS=L1a,L2b,L3,L4,L5,L6&MAP=MFc.map&FORMAT=image%2Fjpeg&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=epsg:29183&WIDTH=612&HEIGHT=516&BBOX=....

The combinatory is:

Many mapfiles (c=1..20, a=1, b=1..5): a templating software generated 20 mapfiles, each with ~5 layers.
One mapfile (c=1, a=1..20, b=1..5): a templating software can generate a mafile with all ~200 layers.


texts: Consuming an asynchronous Geoprocessing service in ESRI's Flex API

texts:
I have been trying to consuming my own asyncronous Geoprocessing service using ESRI's FLEX api, but i am either doing something dumb, or there is some trick involved in getting the results out.

Since it's an asyncronous service, I am using the submitJob()method. The jobComplete event is firing, but the GeoprocessorEvent that i am getting here has a null executeResult.

And if I check using FireBug, I can see that the result is definitely being sent to my browser.

Is there any sample which shows how to get the results  from an asyncronous service?

P.S. I am using ArcGIS server 9.3.1 & Flex API 2.1


texts: MakeTableView of Excel spreadsheet tab gives ERROR 000732?

texts:
I have a very simple model built using Modelbuilder that takes a tab of an Excel spreadsheet and converts it using Make Table view. This works fine in ModelBuilder. I then export the model to Python to give the following code:

import arcpy
Combined_ = "D:\\EAD594\\EUforGIS.xls\\Combined$"
Combined_View = "Combined_View"
arcpy.MakeTableView_management(Combined_, Combined_View,...)

This gives the following error:

  ExecuteError: Failed to execute. Parameters are not valid. ERROR
  000732: Input Table: Dataset D:\EAD594\EUforGIS.xls\Combined$ does not
  exist or is not supported

Any ideas what could be causing this?


texts: Cartodbfying and using ghost tables in CartoDB Enterprise?

texts:
I'm using Carto Enterprise, and I need help on making custom DB tables show up in Carto builder page when you log in.
So far I've turned ON the ghost_tables feature like this.

/opt/carto/builder-tools/embedded/bin/carto-builder-feature-flags.sh enable all ghost_tables

I've also created a table and used cartodbfy.

curl -v -k -X GET https://username.carto.com/user/mbpskt/api/v2/sql?q='CREATE+TABLE+usertable+(+cartodb_id+INTEGER+)'&api_key='key'

SELECT cdb_cartodbfytable('username','username.usertable');

After cartodbfying, it will create the_geom and the_geom_webmercator columns, but it just NEVER shows up in the builder screen.

I think I've done everything possible, but it just doesn't work. Please help me.


texts: Remove/subset raster bands in Python GDAL

texts:
I'm trying to remove bands from a raster I have. Basically I just want the first band as a new GeoTIFF file.

Here is my code so far:

#!/usr/local/bin/python
from osgeo import gdal

# Load my GeoTIFF
tif = gdal.Open("mygeotif.tif")
tif.RasterCount == 3  # should return True

# Create a copy with same projection, etc.
gdal.GetDriverByName("GTiff").CreateCopy("mynewgeotif.tif", tif)
newtif = gdal.Open("mygeotif_band1.tif")

# remove bands 2 and 3 here

# Close files
newtif = None
tif = None

The biggest problem is that I can't find any documentation for writing out only one band or for removing bands.


texts:
I'm looking for a Gazetteer service. I only need the geocoding part and I need two entry points:

name look-up. Takes a query string, check for matching toponyms and
suggest alternatives. Ideally the scope of search could be limited to
certain items, eg only countries, regions and towns. Returns
records with names, boundaries and an id.
get shape. Takes the record id and sends back the record shape coordinates as
WKT (or other). Ideally, it can be simplified before being sent, eg
reduce the number of points.

I have used nominatim until now but they don't allow the get by id, so if I need to download the whole polygon each time, which is less than optimal... also it can only be limited to one feature type at the time.
As an alternative, I could download and install the whole nominatim package and implement the missing features myself.

Before I go into that, I'd like to know if there is anything existing, possibly commercial product, which does what I am looking for?


texts: Where can I download free satellite imagery of Cambodia for offline use?

texts:
I am an entry level GIS user - I have had one semester of GIS in college.

I have successfully downloaded raster satellite imagery from Seamless (tiff files) for my class project.

Now I am trying to develope my GIS skills by volunteering with a non profit group that is doing environmental work with fishing villages in Cambodia.

I would like to download free satellite imagery of central Cambodia - similar to what I did with the Seamless data, to begin my geoprocessing project, any recommendations?


texts:
I'm writing a script to be run on ArcGIS for Desktop, which only runs on Windows machines.

Is there any use for a shebang line at its beginning, and if there is then what should it be (e.g. #!C:\Python27\ArcGIS10.3\python.exe)?


texts: PyQGIS in standalone scripts on Mac, import error

texts:
Following the PyQGIS Developer Cookbook to use PyQGIS in standalone scripts. Fixed the initial errors by setting the PYTHONPATH and DYLD_PYTHON_PATH to the appropriate folders but am now faced with the error:

ImportError: dlopen(/Applications/QGIS3.4.app/Contents/MacOS/../Resources/python/PyQt5/QtCore.so, 2): Library not loaded: @executable_path/../Frameworks/QtCore.framework/QtCore
  Referenced from: /Applications/QGIS3.4.app/Contents/Resources/python/PyQt5/QtCore.so
  Reason: image not found

It loads fine in the QGIS interpreter, and I ran sys.path in QGIS, then ran sys.path.append in my external interpreter for every path that QGIS is connected with from the sys.path results, but no luck. I've found the QtCore.so image but directing the path to that file also doesn't work.

Does anybody know how to get around this issue?


texts: polygons to polygon conversion with attribute name

texts:
1) I have a polygon shapefile with 1000 multi-part polygons in it. Each Row of the table is a multipart polygon.

2) I need 1000 shapefiles with 1 multi-part polygon in each.

3) Each output shapefile with 1 multi-part polygon needs to be named with a name from the attribute table (1) unique to each input polygon.

Open to ArcGIS, ArcPy, Model Builder, OGR, QGIS solutions.


texts: What's the most appropriate way to load Mapbox Studio tiles in Leaflet?

texts:
I'm updating a Leaflet map's basemap, and I'd like to use Mapbox to create this basemap. At first I thought I would have to ditch Leaflet and switch entirely to Mapbox GL. After digging around in the Mapbox API documentation I found an endpoint for requesting raster tiles from Mapbox Studio styles and I got this working in Leaflet.

Is there a problem with loading the tiles this way? Can I expect this method to work for the foreseeable future despite the emphasis Mapbox has been putting on vector tiles? It seems like my only other option is mapbox-gl-leaflet, which is experimental.


texts: Buffers Not Overlapping correctly in ArcGIS for Desktop?

texts:
I have two records in one shapefile that I am trying to create two separate buffers around.

When I go to Geoprocessing - Buffer I choose the desired shapefile, set the desired buffer distance, make the Side Type Full, End Type Round, and Dissolve Type None.

When I run this the first buffer will be an oval as seen in the image, but the second buffer is a weird shape.

How can I make both buffers be an oval?


texts: Calculate NDVI values from images extracted by dates from a list

texts:
Please find attached basic architecture for my data-extraction procedure of NDVI values per pixel, for all images taken in a certain ROI in given dates:

The code below represent the architecture on a small scale of date-range and applied only for SEN2 images (to be done later - the same for SEN3 images):

//Dates of Interest
var start = ee.Date("2018-04-01");
var finish = ee.Date("2018-05-01");

///--------------------- SENTINEL2 Collection ---------------------------------------///
var sen2 = ee.ImageCollection('COPERNICUS/S2_SR')
.filterDate(start, finish)
.filterBounds(region)
.map(function(image){
  return image
  .clipToCollection(region)
  .updateMask(image.select('QA60').not())
})

///--------------------- Unique list creation ---------------------------------------///
sen2 = sen2.map(function(image){
  return image.set('simpleDateMillis', ee.Date(image.date().format('YYYY-MM-dd')).millis())
});
var listMillisSimple = ee.List(sen2.aggregate_array('simpleDateMillis'));
var uniqueDatesSimple = listMillisSimple.distinct();
var SEN2_day0 = uniqueDatesSimple.map(function(date){
  return ee.Date(date);
});

print("SEN2 day0:",SEN2_day0);

//------------------------------ AOI ------------------------------//

//define all AOI
var field1 = ee.FeatureCollection(region.filterMetadata("system:index","equals","00000000000000000000"));

//-----------------------NDVI CALC-----------------------//

//Get the ndvi for the specific field and time
function get_ndvi_sen2(subRegion, image) {
  var ndvi_s2 = image
     .normalizedDifference(['B8', 'B4'])
     .rename('ndvi_s2');

  var ndvi_s2_Features = ee.FeatureCollection(ee.List(ndvi_s2
  .reduceRegion({
    reducer: ee.Reducer.toList(),
    geometry: subRegion,
    scale: 10
  })
  .get('ndvi_s2'))
  .map(function (ndvi_s2) {
    return ee.Feature(null, {ndvi_s2: ndvi_s2});
  })
  );
  return ndvi_s2_Features;
}

//------------------ CSV Export ---------------------//
var sen2_disc, sen2_selc = 'ndvi_s2';
var sen2_folder = 'GOOD/SEN2'

function write_to_drive(col, disc, selc, folder, field){
  print("Col: ",col, "Disc: ",disc, "Selc: ",selc, "Folder: ",folder, "Field: ",field, date);
  Export.table.toDrive({
  collection: col,
  fileNamePrefix: disc+' '+field,
  selectors: [selc],
  folder: folder+'/'+field
  })
}

var feat_per_img = function(sngl_img){
  return get_ndvi_sen2(field1, sngl_img)
}

function write_to_drive_front(feat){
  write_to_drive(feat, 'NDVI_Sen2', sen2_selc, sen2_folder, 'field1');
}

//------ test run for SEN2_day0 list -----//
var allImagesInRange = sen2.filter(ee.Filter.inList("system:time_start", SEN2_day0))

write_to_drive_front(ee.FeatureCollection(allImagesInRange.map(feat_per_img)));

I hope it is not too messy.
As of now - the problem is that I get empty feature collection that needs to be later passed to the function that writes it to CSV file. So I get an empty CSV file.

My main question is - why doesn't the FeatureCollection pass on properly so the Export.table.toDrive could write it to CSV?

Another question is why isn't the code extracting as many CSV files as the number of dates on the list?
 I have mapped the image collection so the procedure of calculating and extracting ndvi values should be on every image of the ImageCollection.

Any ideas?


texts: How to display a .jpg file in ArcMap for georeferencing?

texts:
I'm struggling to georeference a .jpg map. When I load it in to ArcMap, I cannot see the picture at all, so I cannot georeference it. (I tried to zoom to the map layer but the map just does not appear).

Does anyone have a solution for this?


texts: Search by street name/postcode in OpenLayers

texts: Can anyone help me search by street name in OpenLayers? I essentially want to be able to search for streets/postcodes in OpenLayers in the same way I can search on Google maps. Is this possible?


texts: Zonal Statistics as table creates output but some are Null

texts:
I have temperature as a raster file and different ecoregions of a province as polygons.
What I want to do is calculate the mean temperature per ecoregion.

I succeeded in doing this by Zonal Statistics which creates an output raster, but what I actually want, as I realized later, is I the temperature value per ecoregion. So how do I assign that value to the ecoregion?

I thought I can do zonal statistics as table and then join this table to the ecoregion.
It does work but for some reason it does not calculate values for some of the ecoregions (1400 ouf of 7400 get null values). I tried to find a pattern in the ecoregions that didn't get assigned a value, so e.g if they were too small or have a too complicated shape or if they are at the edge of the province...but no..they are just randomly distributed.

Any ideas or something that I am missing?


texts: In the QGIS Browser, is there a way to restrict table deletions/truncations without doing this through SQL Server permissions? As a GIS admin if I want to remove a layer from the map it is too easy to accidentally right click on the feature class in the browser instead, as I have almost deleted a SQL feature class instead of removing the layer. Other GIS users will have this restricted through SQL Permissions.


texts: Oriented minimum bounding box tool missing from QGIS processing toolbox?

texts:
I am working with QGIS 2.14.9 and I have realized that some processes are missing. For example, the Oriented Minimum Bounding Box (OMBB) tool is missing in QGIS vector general tools. It was the same with 2.14.7 version.

I have deleted the processing directory and re-launched QGIS without success.

Is there a way to "update" processing tools?


texts: Determining the cosine of the angles between two lines using dot product

texts:
SCENARIO:I have a certain buffered region (lets say a one way road segment), which will have only one allowed direction represented by vector A. And, there are several other lines which may or may not be in the direction of the road, which are represented by the vector B.

I am trying to find the cosine of the angle between the lines using the cosine(angle)= (A • B) / ||A|| ||B||.

At first, I just wanted to know how to calculate this for one pair of VECTORS A and B. I want to use this to automate the process for the other lines (vector B) in order to find the cosine of the angle.

Is there any tools, so that I can use this ModelBuilder OR is there some code that I can use for the calculation of cosine of the angle?


texts: GPS accuracy deviation the same in an area?

texts: I am using a GPS receiver to perform a speedtest over a few meters. I’ve read that the GPS is accurate up to 4 meters. Is the deviation the same in an open area and also over a short period of time? Or can this vary everywhere and also over a (short) period of time?


texts: Calculate all 'Null' integer fields to '0'?

texts:
Is there an easy way to calculate all Null integer fields to 0 in a GDB?

Selecting Null values and field calculating each column to 0 is not a viable option as I need to do this for a few hundred columns across multiple layers.


texts: What is difference between carto.js and cartoDB?

texts:
I am exploring how to use carto.js map in a web application. I had seen couple of samples where carto.JS is used and some of them used cartoDB by including a supporting js file in the document.

What is the difference between carto.js and cartoDB?


texts: Copy Features tool run after Make XY Event Layer on CSV file fails to complete?

texts:
Using the Make XY Event Layer tool I converted a table to a layer.

I am trying to convert the newly created layer to a feature class in my file geodatabase.

I have been able to do this in ModelBuilder but when I exported the model to a script and ran that it failed to complete - see Export as Python script from ModelBuilder never completes?.

When run from its tool dialog it makes no progress and does not give an error (while in ModelBuilder it takes about 7.5 seconds for it to run).

I tried to convert the layer to a shapefile with ConversionTools>ToShapefile and get the same result.

Finally I right clicked the layer I created with MakeXYEventLayer and performed Date>ExportData>Exported it to my file geodatabase and it worked.

How do I get this to work via the tool dialog (so that it will then work in the Python script)?


texts: st_transform from GRS80 returns empty geometries

texts:
I wish to transform a set of South Korean province shapes from CRS (Coordinate Reference System) GRS80, but st_transform returns empty geometries.

I have downloaded the dataset from this page (in Korean) which specifies the CRS as GRS80. I attach the file, as it is not accessible outside Korea.

Next is the R code resulting in empty geometries:

library(sf)

#Read shapes into R
sca_nsdi <- st_read("Z_NGII_N3A_G0010000.shp", stringsAsFactors=F)

#CRS is not specified
is.na(st_crs(sca_nsdi))

#Assign CRS as EPSG code of GRS80 per this reference: epsg.io/42310
st_crs(sca_nsdi) <- 42310

#All geometries full
sum(!st_is_empty(sca_nsdi))==nrow(sca_nsdi)

#Transform CRS
sca_nsdi <- st_transform(sca_nsdi, 5179)

#Resulting geometries all empty
sum(st_is_empty(sca_nsdi))==nrow(sca_nsdi)


texts: Run Python script with arguments or return values in ArcGIS Workflow

texts:
I found out how to run Python script in custom Step in ArcGIS Workflows. It is possible if I use Executable Execution Type, as executable to run is called python.exe and argument is path to the script. Like this:

But I want to find out how to add some other arguments, like layer names that are created during previous Workflow Steps.
How can I get Return Value from some Step as an argument into this Python script?


texts: I've coded up Karney's transformations from UTM to latitude/longitude coordinates. My work however requires me to know if this conversion is in geocentric or geodetic latitude. I have not been able to find this information online. Does anyone know if Karney's "Transverse Mercator with an accuracy of a few nanometers" is conversions to geocentric or geodetic latitude? I have an educated guess, but could someone point me in the right direction?


texts:
I am working on creating a plugin on QGIS 3.10. This plugin loads two layers in my QGIS project. I would like to zoom in the layer extent of AEE_layer. Unfortunately despite the code that specifies to zoom the canvas on the AEE_layer extent, it does not work. The canvas systematically zooms in on the extent of all the layers that are loaded in my project.
project = QgsProject.instance()
canvas = iface.mapCanvas()

project.write('C:/CartoIni/Projet/Projet.qgs')

groupe_aire = root.addGroup("Aires d'étude")

AEE_layer = QgsVectorLayer("C:/Users/Thysma/Documents/Stage/Data/AEE.shp", 'AEE','ogr')
QgsProject.instance().addMapLayer(AEE_layer, False)
groupe_aire.addLayer(AEE_layer)

reg_layer = QgsVectorLayer("C:/Users/Thysma/Documents/Stage/Data/REGION.shp", 'Limites régionales','ogr')
QgsProject.instance().addMapLayer(reg_layer, False)
groupe_aire.addLayer(reg_layer)

extent = AEE_layer.extent()
canvas.setExtent(extent)
canvas.refresh()


texts:
I am attempting to add a query layer to ArcGIS Pro 2.1.2 via SQL Server.
The table consists of Ticket_ID, Latitude, Longitude. The query layer works if I only add the top 1000 rows. However, when I attempt to add all 79,592 rows, I keep returning the error: Failed to load data. Error: Failed to retrieve a page of rows.

Does ArcGIS Pro have a limit on query layer rows?


texts: How can I use jp2 aerials photos in QGIS?

texts:
I have previously been able to add raster data in the jp2 format to QGIS. I am able to do so on 1.8 installed on my laptop. Unfortunately when I attempt to do this on my desktop, it claims that the file type is unsupported. I installed using the OsGeo installer.

I think I may have missed a setting on the advanced install?


texts: Extracting field value from Info table( .info ) in ArcGIS

texts:
I am working on a Python snippet try to use arcpy.da.SearchCursor() to extract value of certain field from a Info Table which generated by executing arcpy.Statistics_analysis().Here's my code:

arcpy.Statistics_analysis(fc,stats,[["Area","SUM"]],"")
Table = "D:/xxx/stats.info"
with arcpy.da.SearchCursor(Table,["AREA_SUM"]) as cursor:
    for row in cursor:
        Number = row[0]
print Number

However I always got RuntimeError: cannot open 'D:/xxx/stats.info'. I think it's caused by improper file name extention, but I'm not sure what the correct one is.

Plus: is there any other possible way of extracting value from Info Table using arcpy?


texts:
I have a PostGIS layer with linestrings, all linestrings are cut in each line intersection.
The linestrings have no other information than the geometry and an unique id for each linestring. I would like to join linestrings that have a start-/end-point corresponding. I don't care how the lines in eg a 4-way intersection are joined as long as the result is joined linestrings and not multi-linestrings.
Can't find a way to do this.
The layer looks like this:


texts: How to run a python script on the server, returning the results to the client

texts:
Here's the end result I'm trying for:

the user clicks on a feature on an ArcGIS Server JS map
a Python script is triggered on the server
the Python script uses an attribute from the feature (eg its OID) to perform some calculations (using arcpy modules)
the results are returned to the infoWindow as HTML, for display

I guess I'd need to be able to trigger a Python script (containing arcpy commands) on the server, passing parameters between the client and server.

Assuming this is possible, what is the best approach:

via the ArcGIS Server geoprocessing framework
directly via an AJAX call to the Python script (and is that even possible)?
something else?


texts: Offset in pixels observed in outputted raster when "Clip raster by mask layer" of QGIS is used

texts:
I am using QGIS tool "Clip raster by mask layer" in Processing Toolbox for clipping satellite(georeferenced) image raster using shape polygon with "clip to mask layer's extent" setting switched on.
I observed that the output clipped raster's pixels have an offset of around 0.5 pixel size with respect to source image.
Apparently there is a mistake in the algorithm of this tool especially while setting the Geotransform (which decides the positioning) of the output clipped raster.
How can I access the algorithm of this tool?
EDIT:
Found that the pycode ClipByMask.py is apparently (correct me if needed) not the code I am looking for, since it returns argument that calls other GDAL code. Below is the line that shows it:

return ['gdalwarp', GdalUtils.escapeAndJoin(arguments)]

Also observed that the pycode ClipByMask.py has following imports:

from processing.algs.gdal.GdalAlgorithm import GdalAlgorithm
from processing.algs.gdal.GdalUtils import GdalUtils

I was looking for the lines which contain geotransform[0] and geotransform[3] which define the positioning of the image top-left corner, but only found:

      arguments.append(unicode(geoTransform[1]))
      arguments.append(unicode(geoTransform[5]))

above lines are related to only pixel width and height.
I am in search of the GDAL code that has code snippets that involve positioning the raster and has functions to do the raster clipping operation.


texts: distance() function in field calculator gives unexpected results

texts:
I have two layers both containing points - each point in the first layer has a corresponding point in the second layer (same primary key). One layer has CRS EPSG:27700 and the other layer has EPSG:4326. Project CRS is EPSG:27700.

I want to calculate the distance between each pair of points. I tried the Distance Matrix, however the distances returned were nothing like those returned by Measure Line (which I used as a check).

I then tried using the field calculator, with the expression:

distance($geometry,geometry(getFeature('LAYER2','KEY', "KEY" )))

However again the values are much too big, e.g. for one pair of points:

  Measure Line give 99.654m; distance gives 849950.964330714.

I wondered if the problem related to the different coordinate systems, so I saved the EPSG:4326 layer as a shapefile with CRS EPSG:27700. This time distance() returns sensible results, e.g. for same pair of points:

  Measure Line give 99.654m; distance gives 99.6801982600587.

For a given project CRS, shouldn't distance() give the same results irrespective of the coordinate systems of the individual layers?


texts:
In qgis
I have a grid,
 i have a point vector
 for each cell of the grid, i want to make the mean of one attribute of the vector.
How can I do that by not making a raster?


texts: Cropping a Shapefile into a Rectangular Box with Python pyshp

texts:
Using 'pyshp' I am hoping to reduce the set of US highway shapefiles to only those in Puerto Rico.  The idea is very easy since I only need those pieces which fit into the window.

>> import shapefile
>> sf = shapefile.Reader("roadl_usa")
>> len(sf.records())

193504

Not familiar with the writing object but hoping it's easy since I am not building from scratch.  How do I select only those records in the bounding box [-67.5, -65.0] and [17.8, 18.6] and write a new shapefile?

What is the output like?  I am wondering if .shx, .txt, .dbs also get exported and what they mean.

This script almost works except it leads to IndexError: list index out of range - maybe this is because I am deleting things in the process of iterating?

import shapefile
e = shapefile.Editor("roadl_usa")

prbbox = [-67.5, 17.8,-65.0, 18.6]

L = len(e.shapes())
for t in range(L):
    if f( e.shape(t).bbox, prbbox) > 0:
        e.delete(t)

e.save('pr-roads')


texts: Overlay a raster image using Equirectangular projection (lat/lon) onto OSM/Google map with openlayers?

texts:
I have a forecast map of 6x6 lat/lon extend from 104E-110E, 10N-16N as an static PNG image and want to make it an overlay on the basemap of OSM/Google map using openlayers 3. The problem is OSM uses Mercator project while my images uses Equirectangular and the maps don't fit.

You can see my implementation at
http://thoitiet.hus.vnu.edu.vn/?m=meteo&a=ol

Maybe there are methods where I can transform the image to Mercator first, but it's much better if openlayers can handle the correct projection. How/Is it possible that can I do it?


texts: Adjacency matrix for polygons using Java

texts:
I am trying to create an adjacency matrix from a set of polygons. That is, I have a bunch of polygons and I want to identify which polygons have a common edge or "touch" each other. Once I find this information, I want to create an n x n matrix that indicates whether those each polygon either touches or does not touch the other polygon. Something like this:

      Poly1    Poly2   Poly3
Poly1   0        1       0
Poly2   1        0       1
Poly3   0        1       0

That is the goal anyway. I was wondering if anyone knew of a Java package or such that could help me do this. I read some posts on doing this using VBA but I am wondering if Java Eclipse can be used for this purpose because I need to implement the adjacency matrix in a bigger program that I am writing on Java.

Edit:
I have the polygon as a shapefile drawing. I think triangular Matrix will be fine because of the symmetry.


texts: Efficient way for Buffering polygons with a large dataset

texts:
I use artificial grid cells (Polygons) for Sub-Saharan Africa to do a sub-national study. I have many variables assigned to each grid cell. Now, I need to include information (development aid programs) from the surrounding grid cells as control variables. By now, I try it with this code and with the raster package. (gid stands for the cell id)

  for(entries in nrow(dataset)){
      current_cell <- dataset[entries,]
      current_cell_id <- current_cell$gid
      prio_cell <- cellsSSA[cellsSSA$gid == current_cell_id]
      prio_cell_buffer <- buffer(prio_cell,50000)
      neighbor_prio_gids <- cellsSSA$gid[prio_cell_buffer,]
      neighbor_data <- subset(dataset, is.element(gid,neighbor_prio_gids))
      dataset$neighbor_aid[entries] <- min(neighbor_data$AID)
    }

As I want to estimate the effect for many years, I got 85.000 grid cells and R runs now for 6 days. Any ideas to make this code more efficient or any other way to solve my problem? Is there a buffer function that is more efficient?


texts: How to use GDAL in java classes with Eclipse?

texts:
I want to use GDAL in java classes with Eclipse.

I had gdal.jar to the build path but it does not work ("Native library load failed.
java.lang.UnsatisfiedLinkError: no gdaljni in java.library.path").

I try to add libraries (gdalconstjni.dll + gdaljni.dll + ogrjni.dll + osrjni.dll) and to modify environment variables (as some tutorials advocated) but without success.

What should I do to use GDAL in my java classes with Eclipse?


texts: Where can I download Real time Traffic and weather data for the world?

texts: Where can I  download real time data for traffic and weather? I'd like to use it as an overlay on a mapping service of my choice (eg. Google maps, OpenLayers etc.)


texts: How to solve error: global name 'qgis' is not defined

texts:
I use this code as a rubber from my QGIS canvas:

def obrisi_sve(iface):
    rubber = [i for i in iface.mapCanvas().scene().items() if issubclass(type(i), qgis.utils.gui.QgsRubberBand)]
    for ver in rubber:
        if ver in iface.mapCanvas().scene().items():
            iface.mapCanvas().scene().removeItem(ver)

When I run it from Python console inside QGIS it works ok, but when I use it in plugin I get error: global name 'qgis' is not defined. To call this file
from plugin I use these lines of code:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from math import *
import os

import brisanje_rubber

from PyQt4 import QtGui, uic
from PyQt4.QtCore import QObject
from PyQt4.QtCore import SIGNAL
from PyQt4.QtCore import QVariant
from qgis.gui import QgsMessageBar
from qgis.gui import QgsMapToolEmitPoint
from qgis.utils import iface
from qgis import *

def __init__(self, parent=None):
    self.setupUi(self)
    self.pushButton_6.clicked.connect(self.obrisi_IPP)
def obrisi_IPP(self):
    brisanje_rubber.obrisi_sve(iface)

How can I solve this?


texts:
I have a layer with polygons in QGIS. The one selected on the bottom has a hole in it, the hole is also a shape. So I have the donut and his hole. How can I merge these two polygons?

EDIT: I didn't emphasise this but I want to utilise the solution for my whole dataset of around 7 million polygons.


texts: Saving matplotlib map to shapefile in python?

texts:
I already plot a prediction map I made with xgboost modal by matplotlib function and I want to save this map as a shapefile so I can work with it in ArcGIS, how can I do that?

Picture of the map:


texts: How to convert raster (big size) to points?

texts:
I have a raster : population density of Ivory Coast. I need to convert it to points, but I met several problems because of the size I think.

Here is what I tried :

1) GRASS GIS 7 command : r.to.vect => It works but takes several hours !

2) SAGA (2.1.2) : Raster values to points => It doesn't work (Couldn't load gdal_MrSID.dll -> but gdal_MrSID.dll is well installed in the SAGA package).

3) GDAL/ORG -> GDAL Conversion -> gdal2xyz -> It works but takes a lot of time (after 15min Processing algorithm still on 0%).

4) Raster -> Conversion -> Translate -> Save as ASCII grided XYZ. Then Layer -> Add Layer -> Add Delimited Text Layer :

     File name : the created file : xyz file
     File format : Custom delimiters : Tab & Space
     Record option : Unchecked the "First record has field names"
     Geometry definition : Point coordinates : X field = field_1 ; Y field = field_2.

Layer appears in the Layers Panel (contains 40 millions points).

Does anyone have a solution?


texts: How to define QgsVectorLayer by QDomNode?

texts:
I am using QGIS 2.4, Python 2.7.5, win 8.1.

I am trying to read layers from QGIS project. When I try define QgsVectorLayer using QDomNode and readXML function, following error appears:

layer.readXML(layerDOM)
AttributeError: 'QgsVectorLayer' object has no attribute 'readXML'

code is here:

from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import QApplication
from PyQt4.QtXml import *

app = QgsApplication([], True)
print app.initQgis()

doc = QDomDocument()
doc.setContent(QFile("composer.qgs"), False)

layers = doc.elementsByTagName("maplayer")

for x in xrange(layers.count()):
    layerDOM = layers.at(x)
    layer = QgsVectorLayer()
    layer.readXML(layerDOM)

app.exitQgis()

Does anybody know, why does the error appear?


texts: How to use the % Modulo operator to select points for labeling?

texts:
I have a shape file with point numbers for individual ID's. I would like to select point numbers divisible by a number, and if the remainder is 0 then I would like to label that point.

I have tried about 1000 different combinations using the Modulo operator, but I continue to get an "Expression is invalid" error.

Pseudocode for the expression is as follows:

Case
When % ("Point_Number", 7) = 0
Then "Point_Number"
End


texts: Getting Bing imagery to fully load in QGIS?

texts:
I'm using QGIS 2.8 and by using the OpenLayers plugin, I can access Bing aerial imagery.  When I try to load it about half the tiles are missing.

How can I get it to fully load so that my map is completely covered by the imagery?


texts: Extracting shorelines from BW aerial photos using ArcMap

texts:
I'm currently trying to analyze rates of shoreline change for an estuarine reserve by extracting & vectorizing shorelines from historic aerial photos, and then comparing those to the modern shorelines.
All of the literature I've read thus far seems to rely on the use of an IR band, which I don't have access to because a large portion of my data are just BW scans of photos taken by the USDA between 1937-1980. Extracting the modern shorelines isn't an issue because I have access to IR via satellite imagery, but I'm unsure how to approach the years in which that isn't an option.
So far, I've tried Iso Cluster Unsupervised Classification in ArcMap, with mixed results mostly depending on the contrast and quality of the images themselves.
Does anyone here know of a more efficient and reliable way to extract shorelines without an IR band?


texts: How to resolve Projection library error. proj error "Permission denied" in MapServer on IIS

texts:
I'm trying to get MapServer up and running under IIS on my Windows 7 x64 machine.

I've installed the latest 7.0.1 release msi from GISInternals

Added mapserv.exe to the Allowed list in the ISAPI and CGI restrictions

Created a map file, and tested it locally and it generates the expected image

But when I browse to localhost, http://localhost/mapserver/mapserv.exe?map=test.map&mode=map, I get an error ...

msProcessProjection(): Projection library error. proj error "Permission denied" for "init=epsg:4326"


texts: PostGIS | Get degrees clockwise from north of a linestring's beginning and end

texts:
I have a table full of linestrings in PostGIS. I would like to know the degrees clockwise from north of the starting point of the linestring and also of the end point.

I was thinking that maybe creating a circle with both points could help to get this two values (example: X=32º, Y=209º) but i have no clue how to achieve this....

Regards


texts: How should I use GDAL to perform a transformation?

texts:
I am a student who wishes to use GDAL to transform a series of images from Lambert Conformal Conic (EPSG: 9802) to Google Earth's WGS 84 (EPSG: 4326) for a web mapping project.  I know that in order to perform this operation I need to use gdaltransform, which is an executable that can be used via command line.

My issue is not being able to wrap my head around how to use gdaltransform correctly.  I have read several documents that source this link http://gdal.org/gdaltransform.html, but I have been unable to replicate the tutorial results.  Can someone please help me understand what I am missing.

Source Code:
This code is what I typed in replicate the demo found on the above link.

Input:   gdaltransform -s_srs EPSG:9802 -t_srs EPSG:4326

Return:  0 0 0

Input:   177502 311865

Return:  177502 311865 0

Note:  The return should read 244510.77404604 166154.532871342 -1046.79270555763 not 177502 311865 0


texts: Exporting random forest decision tree (or thresholds) Google Earth Engine

texts:
Is there any way to export/view the generated decision tree from a random forest classifier in Google Earth Engine?

I know that

classifier.explain()

works for a CART classifier. Anything similar for

ee.Classifier.randomForest()

Alternatively, is there any way to export/view the thresholds used in the random forest decision tree?


texts: Removing properties from feature in OpenLayers

texts:
I have a feature in my OpenLayers vector layer.
when I do
feature.getProperties()

I get object like this
{geometry: D, drawnAt: "anygeojson.com", date: "", time: "15:20:16", polu: "vgh"}

I want to remove date property from this so that it will look like this
{geometry: D, drawnAt: "anygeojson.com", time: "15:20:16", polu: "vgh"}

How can I do that?


texts:
I have this .asc file of maximum temperature of a local place. This file is a result of vector to raster (GeoTIFF) conversion then GeoTIFF to ASCII grid conversion, both using QGIS. Here is a snippet of the actual file when loaded into WordPad:

ncols        3000
nrows        3000
xllcorner    122.612907409668
yllcorner    9.038180351257
dx           0.000241035461
dy           0.000458419800
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

I need this file as an input for a niche modelling program (DesktopGARP) and when I upload it in the DatasetManager it outputs an error:

> Error loading file contains an invalid ASCII Grid

I'm not sure where the problem is, however I found out that CELLSIZE is needed instead of DX and DY. The first 6 rows, instead, should look something like this:

ncols
nrows
xllcorner
yllcorner
cellsize
NODATA_value

I found a similar question Arc ASCII Grid output problem in QGIS 2.0 but the answer says to use FORCE_CELLSIZE=TRUE. But if I force the DX value into the CELLSIZE value, the raster image would shrink since the values of dx and dy are not nearly the same. Is there a way to manipulate the values without changing the raster's overall size?


texts: Manipulating mpk/gpk/tpk with ArcGIS Runtime WPF

texts:
Is there a way work with those package files within Runtime and NOT in arcmap. It seems very inconvenient for where we have to use these packages, especially for geoprocessing.

Also, is there a way to create/add/remove shape files within those packages?

If anyone has any references, please let me know. This seems to be the trickiest part about Runtime IMO.

Thanks and Regards,
Kev


texts:
Is it possible to monitor phenology of only 2 land cover classes in TIMESAT?
Example is corn and wheat for 7 years. Sowing structure of this crops changes every year, i.e. their location on an image is different every year, different pixels need to be monitored.

I would like to change monitored pixels for every year in order to track the phenology of corn and wheat from 2007 to 2014. Furthermore sowing structure consists of different parcels scattered on an image, approximately 50 MODIS pixels for every year.

I am using MODIS13Q1 product.


texts: Error in displaying WFS layer from MapServer/PostGIS in OpenLayers

texts:
I have read and tested all the relative posts but I still don't have a solution for my problem.

I have a layer suscribers (points) stored in PostGIS database with srid=32232. I would like to display the layer in OpenLayers in geographic projection EPSG:4326. I can see the response with Firebug but get no display.

My code:

map file

MAP
  IMAGETYPE PNG
  EXTENT -180 -90 180 90
  UNITS DD
  SIZE 1500 650
  IMAGECOLOR 255 255 255
  FONTSET "fonts/fonts.list"
  SYMBOLSET "symbols/symbols35.sym"
  WEB
    IMAGEPATH "tmp/tmp_img/"
    IMAGEURL "**/guide/tmp/tmp_img/"
    HEADER "query_header.html"
    FOOTER "query_footer.html"
    METADATA
      "wfs_title" "WFS Guide" ## REQUIRED
     "wfs_onlineresource" "http://*********/cgi-bin/mapserv.exe?" ## Recommended
     "wfs_srs" "EPSG:4326" ## Recommended
     "wfs_abstract" "This is the wfs service of guide" ## Recommended
     "wfs_enable_request" "*"  # necessary
    END
  END
  PROJECTION
    "init=epsg:4326"
  END
  LAYER # suscribers begin here
    PROJECTION
     "init=epsg:32232"
    END
    METADATA
      "wfs_title" "suscribers" ##REQUIRED
      "wfs_srs" "EPSG:4326" ## REQUIRED
      "gml_include_items" "all" ## Optional (serves all attributes for layer)
      "gml_featureid" "gid" ## REQUIRED
      "wfs_enable_request" "*"
    END
    CONNECTIONTYPE postgis
    NAME suscribers
    CONNECTION "user=postgres dbname=postgis20 password=pass"
    DATA "geom FROM public.suscribers USING UNIQUE gid USING SRID=32232"
    STATUS ON
    TYPE POINT
    DUMP TRUE
  END # suscribers ends here
END # end of map

openlayers

var map = new OpenLayers.Map('map',
    { controls: [],
      allOverlays: false,
      tileSize: new OpenLayers.Size(100, 100),
      maxExtent: OpenLayers.Bounds.fromString("-180,-90,180,90"),
      projection: new OpenLayers.Projection("EPSG:4326"),
      displayProjection: new OpenLayers.Projection("EPSG:4326") //mouse position });
var basic = new OpenLayers.Layer.WMS("Basic",
    "http://vmap0.tiles.osgeo.org/wms/vmap0",
    { layers: 'basic' },
    { singleTile: true,
      ratio: 1,
      transitionEffect: 'resize',
      opacity: 0.5
    }
 );
var wfsSuscribers = new OpenLayers.Layer.Vector("States",
    { protocol: new OpenLayers.Protocol.WFS(
        { url: "http://*********/cgi-bin/mapserv.exe?map=C:/ms4w/Apache/htdocs/routing/index.map&",
          featureType: "suscribers",
          featureNS: "http://****/routing"
        }),
    projection: new OpenLayers.Projection("EPSG:4326"),
    strategies: [new OpenLayers.Strategy.BBOX()]
    }
);
map.addLayers([basic, wfsSuscribers]); map.zoomToMaxExtent();
map.setCenter(new OpenLayers.LonLat(11.52, 3.83), 3);

response

<gml:boundedBy>
    <gml:Box srsName="EPSG:4326">
        <gml:coordinates>11.416136,3.714272 11.554812,3.969178</gml:coordinates>
    </gml:Box>
  </gml:boundedBy>
<gml:featureMember>
  <ms:suscribers fid="suscribers.2">
    <gml:boundedBy>
        <gml:Box srsName="EPSG:4326">
            <gml:coordinates>11.521489,3.869029 11.521489,3.869029</gml:coordinates>
        </gml:Box>
    </gml:boundedBy>
    <ms:msGeometry>
    <gml:Point srsName="EPSG:4326">
      <gml:coordinates>11.521489,3.869029</gml:coordinates>
    </gml:Point>
    </ms:msGeometry>
    <ms:gid>2</ms:gid>
    <ms:osm_id>331136682</ms:osm_id>
    <ms:addr_house></ms:addr_house>
    <ms:amenity></ms:amenity>
    <ms:name>Yaoundé</ms:name>
    <ms:barrier></ms:barrier>


texts: Why are the areas covered by 2 Universal Polar Stereographic (UPS) coordinate different?

texts:
According to http://en.wikipedia.org/wiki/Universal_Polar_Stereographic_coordinate_system. the area covering Arctic is less than Antarctic. Any reasons for it?

It is due to the shape of the lands?


texts: Reclassify land use map based on agricultural statistics

texts:
I'd like to create a land cover map of Greece with enhanced information on the spatial distribution of agricultural crops.

As you know, common land cover products (such as ESA CCI LandCover or CORINE) only account for "arable land" without specifying crops cultivated on it. My intention now is to consult agricultural crop statistics of Greek regions in order to create a land cover map including crops. I'd like to distribute the crops over the arable land of each region so that the proportional distribution within each region matches the statistics.

As this is only a statistic approach, the exact location of the different crops remain uncertain, of course. So some kind of random generator distributing a number of, let's say, wheat pixels within a region would be suitable.

Do you have any idea how that could be done using ArcGIS? I have examined the numerous (Re-)Classify and Statistical Analysis tools in ArcGIS, but none of them seem to be suitable. Further I thought about creating a random raster, but there is no way to define custom value proportions, rather some distribution types (normal, gamma, etc.).


texts: Add image not providing north arrow directory in QGIS?

texts:
I'm still using QGIS Brighton to finish a map in Composer and would like to add a north arrow to it.  In previous versions I added the arrow as an image, but when I click on 'add image' and click on the map I only see a cross where the cursor would be. The command history doesn't change and no directory appears from which I can choose and customise an image.

How do I access the image directory?


texts: Populate field with text from another field based on value

texts:
I had two layers (datasets). I have run the tabulate intersection to get a table of intersected features.

This is the table

Now I want to populate field with values based on "BROJ_UP" field, like this, and do that for every field. For example, value from "BROJ_UP" is 1.5, and values from NAM1 are S, Z, Z1, etc. I want to populate field NAM with text from NAM1 field, but based on value from "BROJ_UP" field.

Even better, it could be great to delete duplicate fields, so I after calculations I want to have just one field. And for every different value in "BROJ_UP" field.


texts: Distance from centre point to other points using QGIS?

texts: I have a question of something I actually thought was pretty simple...I have a few points around a centre point. Now I want to know the distance from each point to the centre point. Is it possible to do this without using the simple measurement tool. I mean, is there something that gives me "automatically" all the distances with just one query.


texts:
I try to load a CSV file into QGIS. The file contains a list of points in decimal degree format without E, W, S, or N letters, like 35.876528, 045.982671.

When I add it as a delimited layer, it loads all the points, but on a very small scale near point zero south of Ivory Coast. The distribution of the points are correct on that scale, but are too small to be correctly aligned with the underlying map.

I use the same CRS as the project's. WGS84.

image of the text


texts: How to extract data from a POINT column in PostGIS using JAVA (jdbc)

texts:
How do you get a value from a column which has a point data type.

I've just switched over from MySQL to PostGIS. I have a column which is called point and is of type Point.

I want to get that value as a result of a SQL query.

I've done things like rs.getString("text") and rs.getInt("number") for other columns and they work fine, but doing

rs.getObject("point")

Causes me to get the error

org.postgresql.util.PSQLException: This connection has been closed.

Here's more code, I thought it might a real-dummies question so I didn't stick it in before.

Essentially I have a value in my application which can take on a variety of data types. Not having a truly OO database, I've made a table with lots of columns, depending on which is is populated, I will know what datatype the value is for that record.

public static Set<STriple> readStatements(String subUri, String preUri) {
    Set<STriple> statements = new HashSet<STriple>();
    String sql = "SELECT sub_id, sub_uri, sub_value, pre_id, pre_uri, pre_value, obj_id, obj_uri, obj_value, number, text, file, url, point FROM triples WHERE sub_uri='"
            + subUri + "' AND pre_uri='" + preUri + "'";
    ResultSet rs = Dao.executeQuery(sql);
    try {
        while (rs.next()) {
            STriple s = new STriple();
            s.setSubject(new SResource(rs.getString("sub_uri")));
            s.setPredicate(new SResource(rs.getString("pre_uri")));
            if (rs.getObject("obj_id") != null) {
                s.setObject(new SResource(rs.getString("obj_uri")));
            } else if (rs.getObject("point").toString() != null) {
                Object pointObject =  rs.getObject("point"));
            } else if (rs.getObject("number") != null) {
                s.setObject(new SNumber(rs.getDouble("number")));
            } else if (rs.getObject("text") != null) {
                s.setObject(new SText(rs.getString("text")));
            } else if (rs.getObject("url") != null) {
                s.setObject(new SURL(rs.getString("url")));
            }
            statements.add(s);
        }
    } catch (SQLException e) {
        e.printStackTrace();
    }
    return statements;
}