texts: Empty TIFF files while exporting images using Google Earth Engine

texts:
I am trying to export aerial images of a city in Japan. I have specified a polygon geometry of the location and also specified it as region in export. But I am getting blank TIFFs in the drive as output. What could have gone wrong?

print('started running script');
var landsat8_col = {
  l8_reflect_tier1:ee.ImageCollection("LANDSAT/LC08/C01/T1_SR"),
  l8_reflect_tier2:ee.ImageCollection("LANDSAT/LC08/C01/T2_SR"),
  l8_atmosphere_tier1:ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA"),
  l8_atmosphere_tier2:ee.ImageCollection("LANDSAT/LC08/C01/T2_TOA"),
  l8_raw_tier1:ee.ImageCollection("LANDSAT/LC08/C01/T1"),
  l8_raw_tier2:ee.ImageCollection("LANDSAT/LC08/C01/T2")
}

var landsat7_col = {
  l7_reflect_tier1:ee.ImageCollection("LANDSAT/LC07/C01/T1_SR"),
  l7_reflect_tier2:ee.ImageCollection("LANDSAT/LC07/C01/T2_SR"),
  l7_atmosphere_tier1:ee.ImageCollection("LANDSAT/LC07/C01/T1_TOA"),
  l7_atmosphere_tier2:ee.ImageCollection("LANDSAT/LC07/C01/T2_TOA"),
  l7_raw_tier1:ee.ImageCollection("LANDSAT/LC07/C01/T1"),
  l7_raw_tier2:ee.ImageCollection("LANDSAT/LC07/C01/T2")
}

var l8_dates = [['2015-01-01','2018-05-31']]

var l7_dates = [['2008-01-01','2008-12-31']]

var point_of_interest = ee.Geometry.Point([137.95,35.72])

var polygon_of_interest = ee.Geometry.Polygon([
          [
            [
              137.93982982635498,
              35.728313419207375
            ],
            [
              137.94046819210052,
              35.728313419207375
            ],
            [
              137.94046819210052,
              35.729663398270645
            ],
            [
              137.93982982635498,
              35.729663398270645
            ],
            [
              137.93982982635498,
              35.728313419207375
            ]
          ]
        ]);

var l8_keys = Object.keys(landsat8_col)
var l7_keys = Object.keys(landsat7_col)

var l8_values = Object.keys(landsat8_col).map(function(key) {
    return landsat8_col[key];
});
var l7_values = Object.keys(landsat7_col).map(function(key) {
    return landsat7_col[key];
});

function getFilteredImages(img_collections,dates,point_of_interest){
  var filtered_imgs = [];
  img_collections.forEach(function(img_collection){
    var img = img_collection.filterBounds(geometry).filterDate(dates[0][0],dates[0][1]).sort('CLOUD_COVER').first();
    filtered_imgs.push(img);
    print(img);
  });
  return filtered_imgs;
}

function export_images(imgs_col){
  for(var i=0;i<imgs_col.length;i++){
    Export.image.toDrive({
      image: imgs_col[i].select(['B4','B3','B2']),
      description: l8_keys[i],
      scale: 1000,
      maxPixels:1e13,
      region:polygon_of_interest,
      fileFormat: 'GeoTIFF',
      folder:"2"
    });
  }
}

var filtered_imgs = getFilteredImages(l8_values,l8_dates,point_of_interest);
var vizParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max:1000,
  gamma: [0.95, 1.5, 2]
};

export_images(filtered_imgs);
print('Done');


texts: St_intersects returning false results in postgis

texts:
I want to intersect point geometry with raster. In postgis I'm executing

  select st_intersects(geom,rast)
  from testraster1,pft_pft

and the result is false even though the points are sitting on raster and both have same SRID i.e. 4326. What is the problem?


texts: How to design labels which look like a fraction (with a fraction line)?

texts:
I want to create labels like a fraction using data from two columns. I have tried to use:

"Field" || '\n' || '――' || '\n' || "Field2"

but it doesn't look good. How can I do it using other method?

Here I have added the picture


texts:
I have a circle and a polygon that I'd like to move to the user's local position. I managed to move the circle latLng but I can't find a way to do the same thing with a polygon.

Can someone help me please?


texts: Problems with matching layers with different projections

texts:
I have two layers with the following projections:

Monte Mario IT1
UTM*

Anyhow, I followed the instructions (open the project, set the projection, select on the fly projection, load the layers but no way, the two layers does not match each other.

I could only obtain the result of picture 1 where in red are the transport flows, in blue the transport network.

The crs of the two layers are:

  red--> +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0
  +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +units=m +no_defs blue--> +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=500053 +y_0=-3999820 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs

If I set one of the layer with the same CRS as the other, the two networks becomes very far from each others (see picture 2)

Which do you think might be the problem?


texts: Qgis & processing doesn't links with Grass, r or saga

texts:
After having formatting my computer (W7-64bits), I have re-installed QGIS 2.0.1-64 bits on it .
But I can not use GRASS, SAGA or R via Toolbox (or Processing). I have checked the path in the Options, and I have tried to uninstall and re-install it again, but it doesn’t work yet.
In the beginning of the re-install, I have a message telling that I QGIS is already installed, but it didn't remain any folder or files link with Qgis.
I checked as well in the option and plugin setting, but everything seem to be fine (plugin activate, plugin repositories connected...).

Did I install something else than QGIS? did I miss a step, etc.?

Many thanks for help!


texts:
I am using Leaflet library to render states shapes. I am using OSM for background. I am setting noWrap:true for tile layer as I want to see only one background map covering entire world. Now I would like to make sure my background map is always fitting to the map div, without showing empty space.

Initial it fits to the layer bounds. But once user pans or zooms, we can see white space as the background image is no more available. If the background image is not available, I want to prevent, pan or zoom. How can I fit the background image completely to the div.

Here is my code:

`


texts: Seeking full list of ISO ALPHA-2 and ISO ALPHA-3 country codes?

texts:
I'm searching for a full list of countries with their country codes.

Like on this page (need full and valid):

http://www.nationsonline.org/oneworld/country_code_list.htm


texts: How do I deal with multiple kml files for one map?

texts:
I'm looking for general guidance on how to approach a simple map.
I have a spreadsheet of seven schools, with fields such as "number of students" and "school budget".
The catchment area for each school is its own KML file.
What I want to do is build a map showing the seven catchment areas in my city, with pop-ups showing the school data for each area.
What's the best way to structure this?


texts:
I'm trying to use gdal/ogr to join a shapefile called myshape.shp with a csv called mycsv.csv.
The join field in the shape in joinshp and the join field in the csv is joincsv.

I've found this command on other post(link) :

ogr2ogr -sql "select inshape.*, joincsv.* from inshape left join 'joincsv.csv'.joincsv on inshape.GISJOIN = joincsv.GISJOIN" shape_join.shp inshape.shp

I tried this but I can not get it to work:

ogr2ogr -sql "select myshape.*, mycsv.* from myshape left join 'mycsv.csv'.mycsv on myshape.joinshp = mycsv.joincsv" shape_join.shp myshape.shp

Can someone help me write this command for my example ?


texts:
I am not getting map legend in GWC WMTS or WMS service, QGIS returns an error "Returned legend image is flawed". but when I am accessing that layer using WMS in QGIS, Map legend is showing with the layer. Also, When I am accessing a small area with the same style, Map legend is showing in GWC WMS service.
I am using GeoServer version 2.12.2.
GWC WMS link: http://localhost/geoserver/gwc/service/wms?
WMS link:  http://localhost/geoserver/wms?

Anyone knows why Map legend is not working in GWC WMS and WMTS?


texts:
Is there anyway to alter elevation of DEM where there are bridges and dams for the flood routing of the River?

I have to calculate discharge of the River from ancient times so at that time there were no dams and bridges.


texts:
I want to import a geonames file. I downloaded the geonames .txt file from http://download.geonames.org/export/dump/

I used:
v.in.geonames input=PH.txt output=ph

but it displayed this error:
ERROR: File 'C:UsersRhenDesktopPHPH.txt' not found

I don't know why it cannot be imported..


texts: Loading specific area from an osm.pbf file into Maperitive

texts: Is there a way to download a specific longtitude / latitude boundary area (a city) from a country osm.pbf file in Maperitive? In my case I am trying to load Greater Sydney from http://download.geofabrik.de/australia-oceania/australia.html I would like to use Maperitive because I have a particular .mrules I like to use


texts:
Under QGIS 1.8 on Windows 7 I have this problem with my pyArchInit plugin: I need to load on the map canvas my geometries from a PostGIS layer. Till the 1.7.4 version I use the method explained on the tutorial.
But I have to load also some Potgres view, and QGIS 1.8 need to know the column that I would like to use as ID.
How can I pass this value with PyQGIS and restore my plugin?
uri = QgsDataSourceURI()
# set host name, port, database name, username and password
uri.setConnection("localhost", "5432", "dbname", "johny", "xxx")
# set database schema, table name, geometry column and optionaly subset (WHERE clause)
uri.setDataSource("public", "roads", "the_geom", "cityid = 2643")

vlayer = QgsVectorLayer(uri.uri(), "layer_name_you_like", "postgres")


texts:
For some unknow reason, my average of 3 hillshades is very dark :

gdaldem hillshade input.tif hillshades_C.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_B.tmp.tif -s 111120 -z 5 -az 355 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_A.tmp.tif -s 111120 -z 5 -az 275 -alt 60 -compute_edges
gdal_calc.py -A hillshades_A.tmp.tif  -B hillshades_B.tmp.tif -C hillshades_C.tmp.tif --outfile=./hillshades.tmp.tif --calc="(A+B+C)/3"

I did tried --calc="(A+B+C)", it get lighter, but still too dark and not as expected.


texts: Non-Image Polygon fill (e.g., diagonals) in Mapnik?

texts:
Some of the output produced by our team  - primarily on updates for 'legacy' projects - is styled with diagonal line fills (an example appears below).

In mapnik, this can be achieved using

<PolygonPatternSymbolizer file="[appropriate image]"/>

The exact same effect is possible "image-free" using background-image in CSS3, styled with appropriate gradients.

My preference is for the latter.

My searches indicate that this is not possible in mapnik i.e. it appears that mapnik's PolygonSymbolizer has parameters (fill, fill-opacity, gamma); PolygonPatternSymbolizer takes (file).

I do notice that PolygonSymbolizer.fill 'accepts expressions', as does PolygonPatternSymbolizer.file - so the question is... can expressions be used to achieve the same effect as the image above?


texts: Converting radians to arithmetic degrees

texts:
I need to convert angles that are in radians from Due east to arithmetic degrees for
Creating arrows based on GPS velocities to show displacement

The Radians definition for my data is:

Degrees definition:

I tried using "degrees = !radians! * 180 / math.pi" in field calculalotr

But this gives negative values which are no higher than -/+ 180


texts: Table Relate Spatial and Tabular data in QGIS 2.0

texts:
I am self taught in GIS, so I apologize ahead of time for my noobness.

In ArcGIS you have two options to associate data in a table to spatial data. Based upon my understanding, the first is Join which allows the join of 1:1 or many:1. Then there is Relate which allows 1:many.

I downloaded soils data from Web Soils Survey. In this data there is a spatial data, a shape file named soilmu_a_ny071.shp. And there is tabular data named component in soildb_US_2003.mdb.

I need to relate the spatial data to many attributes in the tabular data based upon the common field MUKEY. Is there a way to Relate in QGIS 2.0? I also have Postgresql and PostGIS installed, so is there a way to do it through there?


texts: Creating ellipse with SatScan output using QGIS

texts:
I need to use QGIS to create an ellipse employing a SatScan output with information on spatial risk clusters for the occurrence of an.

The output contains the coordinates of the centroid point of the ellipse, as well as E-MINOR, E-MAJOR, E-ANGLE and E-SHAPE.

I've already created a circular buffer using MMQGIS plugin, but I could not do the same for an elliptical shape. Any ideas?


texts: Is it possible to label a postgis Query Layer in ArcMap?

texts:
I've created a PostGIS query layer and I can't label it. I'm trying to work out if I'm dealing with a bug that I can't workaround, a bug that I can workaround, or if it is operator error.

After getting my table and columns to lower case and making sure my data types are supported, I was able to add a query, validate, pick the OID column and add it to ArcMap. I can view the attribute table and identify attributes (point) with no problem.

ArcMap crashes (no error--just disappears) when I turn on labels, which are coming from  a character varying(400) field.

Has anyone run into this and knows what is going on?

Here's what I'm running:

Client:
ArcMap 10.2.1 10.2.1.3497

Database:
"PostgreSQL 9.3.5, compiled by Visual C++ build 1600, 64-bit  ,  PostGIS:2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"


texts:
I've forgotten the name of a particular GIS principle and need a refresher...

Scenario: Let's assume there is some polygon which represents some administrative area. Let's say it's a square, for the sake of simplicity. Let's also assume there are 9 houses somewhere within this square, but we don't know where they are. If I divide this square into 9 equal parts (like a tic-tac-toe board), I want to calculate the number of houses within each square.

Now, clearly, nine houses divided into nine small squares is one house per square. However, I know that all nine houses could be anywhere. All nine of them could be in the top right square like this:

Or there could be three each in the bottom row and none in the top six like this:

If I were to assume one house per square like this,

what is this incorrect assumption called?  Is this the Modifiable Areal Unit Problem (MAUP)?


texts: Mapserver NODATA affects valid raster data pixels (zero)

texts:
I have NOAA nightlight rasters cropped to country boundaries that I am trying to display using Mapserver. As such, I am trying to have NODATA pixels i.e. those outside the country boundaries be transparent. The only way I have got this to work is by combining:

OFFSITE 0 0 0

with

PROCESSING "SCALE=AUTO".

This works well in that the pixels outside the countries do not appear. However, the raster data pixels now change transparency "randomly" as you zoom in. I am guessing this is the SCALE=AUTO at work - it is possible to have valid zero data values.

I should mention that the NODATA reported by gdal is -3.4e+38 and Gimp shows the pixels outside the country boundaries have RGB values -500.0 -500.0 -500.0.

I have tried setting OFFSITE -500 -500 -500. This does not work.

I have no clue how to use -3.4e+38 in the mapfile for NODATA.

Below are a few screenshots as I zoom in showing the shifting "optics".

At default view all seems okay ...

First zoom in shift occurs. Maybe zero pixels in raster?

Zoom in again seems smaller patch

Zoom in again. Different shift

My mapfile is as below:

  MAP
    IMAGETYPE   PNG
    EXTENT          -180 -90 180 90
    SIZE                    800 600
    IMAGECOLOR  "#ffffffff"
    TRANSPARENT ON
    SHAPEPATH     /path/to/nightlights/rasters

    PROJECTION
      "init=epsg:4326"
    END

    OUTPUTFORMAT
      NAME      GEOTIFF
      DRIVER    GDAL/GTiff
      MIMETYPE  image/tiff
      IMAGEMODE PC256
      EXTENSION tif
    END

    WEB
      METADATA
        "wms_title" "Nightlight Rasters"
        "wms_onlineresource" "http://localhost/cgi-bin/mapserv?map=test.map"
        "wms_description" "nightlights"
        "wms_name" "Nightlights"
        "wms_label" "Nightlights"
        "wms_srs" "EPSG:3857"
        "wms_extent" "-180 -90 180 90"
        "wms_formats" "GEOTIFF"
        "wms_enable_request" "*"
      END
    END

    LAYER
      NAME NL_CTRYCODE_VIIRS.M_NLPERIOD_VCMCFG-MTSALL-MEAN-RGFT_GADM-3.6-SHPZIP
      STATUS ON

      TILEINDEX NL_CTRYCODE_VIIRS.M_NLPERIOD_VCMCFG-MTSALL-MEAN-RGFT_GADM-3.6-SHPZIP_tindex.shp
      TILEITEM  "location"
      OFFSITE 0 0 0

      TYPE RASTER
      PROCESSING "SCALE=AUTO"

      PROJECTION
        "init=epsg:4326"
      END

      METADATA
        "DESCRIPTION" "nightlights"
        "wms_title" "NL_CTRYCODE_VIIRS.M_NLPERIOD_VCMCFG-MTSALL-MEAN-RGFT_GADM-3.6-SHPZIP"
        "wms_timeitem"  "nlPeriod"
        "wms_timeextent"  "201401/201412"
      END

      # CLASSITEM "[pixel]"
      #   CLASS
      #     NAME "NODATA"
      #     EXPRESSION ([pixel] = -339999999999999996123846586046231871488)
      #     #EXPRESSION ([pixel] < 0)
      #     STYLE
      #       OPACITY 0
      #     END
      #   END #END CLASS
      # END #END CLASSITEM

    END #END LAYER

  END #END MAP

I have tried using classes but not sure how to do this. Plus it seems pretty expensive in terms of processing. Though I may need to look at this in future to have the bright spots stand out against the dark areas. I can supply a sample raster as well if it will help.


texts: Preventing polygon symbology from rotating with data frame when switching between map sheets in map series?

texts:
I am using ArcGis 9.3, Windows XP and want to know if it's possible to have symbology in my polygons not rotate with data frame when I switch between mapsheets in my map series?

I do NOT have MPSAtlas or PLTS. I DO have DS Mapbook installed.


texts: Adding attribute field to raster layer in ArcGIS Desktop?

texts:
I am from Vietnam and I am still a student. I have to study more about ArcGIS.

How do I add an attribute field in this case?

My project is using Landsat 8 picture to establish a drought map. After calculating LST and VHI I reclassify it and my teacher tells me to add an attribute field to it. I don't know how to do it, I watch people do it in YouTube but when they open the attribute table it is very different.

How can I add an attribute field to my project?

I have to calculate the area of each region too.

How do I do that?


texts: Convert XY (SRID 4326) to UTM in spatialite returns NULL

texts:
How can I convert EPSG4326 point to UTM in spatialite?

SRID=4326;POINT(51.2506238 35.70442974)

Executing X(TRANSFORM(Geometry, 32639)) as X  always returns NULL. What's the proper way of extracting northing and easting from geometry type 4326?


texts: Creating a dynamic KML using a websites API with Python

texts:
I'm trying to create a KML that will dynamically update using an API call in python based off where the user is currently looking. I've tried to include as much commented work as possible for what I've already done, hopefully it's not too much.
Right now, I'm using a given list of locations and the Wikipedia API to (1) Look up the coordinates for the given locations, (2) Parse the JSON output for the Coordinates and the PageID, and (3) Creating a KML that shows the Location and the Webpage linked to the location. This all looks like the following, just for reference:
import requests # Used to access webpages
import re # Used to select and clean text
import objectpath # Used to parse JSON response cleanly
import wikipedia # Used to interact with Wikipedia
import simplekml # Used to create the kml from the csv

# Create url parameters for Wiki API; location names goes between these two
base_URL = "https://en.wikipedia.org/w/api.php?action=query&prop=coordinates&titles="
end_URL = '&format=json'

# Mock list of locations:
locations = ['Colorado', 'Texas', 'Florida']
locations_Meta = []

# Create list of lists that will include the metadata for each location
for x in locations:
    locations_Meta.append([x])

# Get the JSON URL from each location
location_API_Pages = []

for place in locations:
    API_URL = base_URL + str(place) + end_URL # API page to get the coords
    location_API_Pages.append(API_URL) # Add API URLs to list

# Get the coordinates from each locations JSON output
tempList = []
tempList1 = []
url_pre = 'https://en.wikipedia.org/?curid=' # Adding a pageID to this makes a complete URL
count = len(locations) # Length of locations so script knows how many times to run the loop

# Generate the JSON output from the webpage and pull out the coordinate information
for x in location_API_Pages:
    r = requests.get(str(x)) # Access the Wikipedia API
    data = r.json() # Turn response into JSON
    tree_obj = objectpath.Tree(data) # Turn the JSON output into a objectpath tree
    tempList.append(list(tree_obj.execute('$..coordinates'))) # Search the tree for specifically
    tempList1.append(list(tree_obj.execute('$..pages'))) # Search the tree for the pagesIDs

# Seperate the coordinates from the rest of the JSON file
for x in range(0, count):
    for k,v in tempList[x][0].items(): # The piece holding the coords is in a dict
        if k in ('lat', 'lon'):
            locations_Meta[x].append(v) # Append the coords

# Seperate the pageID from the rest of the JSON file and
for x in range(0, count):
    for k,v in tempList1[x][0].items(): # The piece holding the pageIDs is in a dict
        locations_Meta[x].append(url_pre + k) # Append pageIDs to later turn into URLs

# Create KML from the CSV that was created as proof of concept
kml = simplekml.Kml()

for x in range(0, count):
        kml.newpoint(name=locations_Meta[x][0], description=locations_Meta[x][3], coords=[(locations_Meta[x][2],locations_Meta[x][1])])
kml.save('Locations.kml')

(This is my first time parsing JSON with python - or anything. I'm sure there's better ways to do my above code, but that's not what I'm here to ask.)
Now, I know you can dynamically request pages with an API call using a lat and long. What I would like to do is use Wikipedia's built-in search method to dynamically retrieve and display pages from their API call found here:
current_location_lat = 'Somehow get users current view from the loaded KML platform'.
current_location_long = 'Somehow get users current view from the loaded KML platform'

Wiki_API_URL = "https://en.wikipedia.org/w/api.php?action=query&list=geosearch&gscoord={0}%7C-{1}&gsradius=10000&gslimit=100".format(current_location_lat, current_location_long)

I can create a static KML from a fixed data source, but I have no idea how to create this kind of dynamic 'linked' KML. Can anyone help point me in the right direction?


texts: How do I make this subquery factoring work with st_geometry predicates?

texts:
I have the following SQL:

WITH w AS (select shape from schema.feature_class_polygon fcp where fcp.field = 123456)
SELECT count(*)
FROM schema.feature_class_lines fcl
WHERE sde.st_intersects(w.shape, fcl.shape) = 1;

When I run the above SQL, I get the following error:

ORA-00904: "W"."SHAPE": invalid identifier
00904. 00000 -  "%s: invalid identifier"
*Cause:
*Action:
Error at Line: XX Column: YY

My question, how do I make the subquery factoring work with the st_geometry predicates?

The use of subquery factoring works in this ST_Geometry example, albeit not part of a WHERE clause:

WITH tmp AS (select sde.st_centroid(w.shape) shape from schema.feature_class_polygon w where w.field = 123456 )
SELECT sde.st_x(tmp.shape), sde.st_y(tmp.shape) from tmp;


texts: I have been researching some Tennessee Valley Authority (TVA) documents from 1947 that show the location of an item of interest to me as X=3,160,900 and Y=730,000. But, when I try to get Google Maps to display this location, it doesn't like my coordinates. Perhaps TVA's X and Y coordinates were created using an antiquated values system. Can you explain how I can convert TVA's coordinates to modern GPS values that Google Maps will accept?


texts: How do I convert TIFF files to BIL files with header?

texts: In my project, I have hundreds of GeoTIFF files that needed to be converted to .bil files with .hdr (header). They are satellite rainfall data. Can anyone help with a quick method?


texts: Are there any open source tools for the processing of UAV images?

texts: I have aerial photographs taken by a Unmanned Aerial Vehicle (UAV). I want to do ortho mosaicing and DEM creation to do volumetric analysis. Are there any open source tools that process the UAV images?


texts: Creating outline shapefile georeferencing existing map/pictures

texts:
I need to create shapefiles (that need to be adapted and transfered for other programs) with the district outline of a couple of Caribbean countries.

Is georeferencing (using the Georeferencing toolbar) of pictures/maps online a good solution?

Is there any other method I can approach, possibly more professional?

What other measure should I take in order to make the shapefiles created can be used after (without basemap)? (Only thing off is Coordinate system, which in my case is WGS 1984).


texts: Get ST_SummaryStats from Raster Table using geoJSON polygon

texts:
I have a raster dataset imported into POSTGIS with SRS 3857. This table is named as prcp0000_normal3857.

The dataset represents statewide precipitation in Montana and can be downloaded below.
ftp://mco.cfc.umt.edu/prcp/Daymet/yearly_normals/geotiff/PRCP0000_normal.zip

After downloading the data was reprojected to 3857, tiled and renamed prcp0000_normal3857 using the GDAL command below.

gdalwarp -co "TILED=yes" -of GTiff -t_srs epsg:3857 -r bilinear prcp0000_normal.tif prcp0000_normal3857.tif

The data was loaded into the POSTGIS database test1 using the command below

raster2pgsql -I -C -s 3857 /dataFolder/prcp0000_normal3857.tif -F -t 100x100 | psql -d test1

I also have a geojson polygon in SRS 4326 that I would like to use to get some zonal statistics on the raster

The geojson looks like this

'{"type":"Polygon","coordinates":
             [[[-108.797607421875,47.73193447949174],
             [-109.31396484375,47.234489635299184],
             [-107.9296875,46.800059446787316],
             [-106.89697265625,47.12247581664114],
             [-107.02880859375,47.76148371616669],
             [-108.797607421875,47.73193447949174]]]

}'

Based on suggestions from John's answer I am trying the POSTGIS command below

WITH clip(geom) AS
  (SELECT
     ST_Transform(
       ST_SetSRID(
         ST_GeomFromGeoJSON(
       '{"type":"Polygon","coordinates":
           [[[-108.797607421875,47.73193447949174],
                 [-109.31396484375,47.234489635299184],
                 [-107.9296875,46.800059446787316],
                 [-106.89697265625,47.12247581664114],
                 [-107.02880859375,47.76148371616669],
                 [-108.797607421875,47.73193447949174]]]
             }'),
        4326),
     3857)
   )
 SELECT ST_Clip(rast, geom, true)
   FROM prcp0000_normal3857, clip;

POSTGIS just returns the errors below with no other info

I feel like I am close but not indicating the band and perhaps messing up the ST_Transform command

Update: I can now successfully run the ST_Summary stats on the table and make the clip. I just need to know how to bring two together.

The ST_SummaryStats quickly returns the mean on the entire table by using the comment below. I just need to know how to pass the clipped raster piece instead of the entire "prcp0000_normal3857" table

SELECT (ST_SummaryStats(rast)).mean FROM "prcp0000_normal3857";

Some more updates. I've fixed projection issues etc and can successfully run ST_SummaryStats using the intersect between the raster and the geojson. That looks like this

  SELECT (ST_SummaryStats(rast)).*
  FROM "prcp0000_normal3857"
  WHERE st_intersects(rast, (ST_Transform(
    ST_SetSRID(
      ST_GeomFromGeoJSON(
        '{"type":"Polygon","coordinates":
            [[[-116.8402862548828,37.59464778787345],[-116.75582885742188,37.59464778787345],[-116.75582885742188,37.67050562029702],[-116.8402862548828,37.67050562029702],[-116.8402862548828,37.59464778787345]]]
          }'),
     4326),
  3857)));

This returns summary stats below but not exactly clipped since I used st_intersects. Results are below. I am thinking there are multiple rows because my raster is tiled.

Bringing it all back together. When I try and clip the raster and run st_summarystats, I still get an error that the raster and polygon do not intersect. At this point, the data are correctly projected and I feel I am missing something in the syntax of this command.

WITH clip(geom) AS
  (SELECT
     (ST_Transform(
       ST_SetSRID(
         ST_GeomFromGeoJSON(
           '{"type":"Polygon","coordinates":
               [[[-116.8402862548828,37.59464778787345],[-116.75582885742188,37.59464778787345],[-116.75582885742188,37.67050562029702],[-116.8402862548828,37.67050562029702],[-116.8402862548828,37.59464778787345]]]
             }'),
        4326),
     3857))
   )
SELECT (ST_SummaryStats(ST_Clip(rast, geom, true))).* FROM
"prcp0000_normal3857", clip;


texts: How to map our data with Google Map data

texts:
I want to add all the bus route data and create a route finding app for my country since Google map has no data about bus routes in my country but with JavaScript API i could get the shortest path between two cities. e.g if i am going to D from A. it says from A to B then turn left go to C and then go to D likewise.

My idea is to get those points and i have my database with bus stops and route numbers, map the data and create my app. but the problem is the points which gives by Google map and the bus stops are not always the same. so i want to know whether there is a way in Google maps to add my bus stop data and bus route information along a particular road then i think it can easily be searched.

Thanks..


texts: GeoJSON and OSM layers with different projection

texts:
I have a GeoJSON from NaturalEarth for world state boundaries and I want to display it with only without backgroung. But also, I want to add a OSM layer as another baselayer option.

The code I am using is the following one, but when I check it, as OSM is in another projection, the center map does not work. However if I change the projection of "mapDiv" the one that does not work is the vector layer with the GeoJSON information.

    <script type="text/javascript">
    var lon = 5;
    var lat = 40;
    var zoom = 5;
    var map;

    function Initialize(){
        // Map
     map = new OpenLayers.Map('mapDiv', {
            controls:[
                new OpenLayers.Control.Navigation(),
                new OpenLayers.Control.PanZoomBar(),
                new OpenLayers.Control.LayerSwitcher()]
            } );

        // OSM
        var osmLayer = new OpenLayers.Layer.OSM("OSM");
        map.addLayer(osmLayer);

        // GeoJSON
        var geojsonLayer = new OpenLayers.Layer.Vector("GeoJSON", {
            isBaseLayer: true,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "world_0.json",
                format: new OpenLayers.Format.GeoJSON()
            })
        });
        map.addLayer(geojsonLayer);

        // Options
        map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);
    }

</script>

I have try something similar to this example, but I continue having problems when switching between the layers.

The GeoJSON has been created with SHP files from NaturalEarth with ogr2ogr command.


texts:
I am trying to create a google earth like wep application. My problem is that I can't get to render a full view of the earth.

The zoom_all() function generates the following error :

RuntimeError: could not zoom to combined layer extents using zoom_all because proj4 could not back project any layer extents into the map srs (set map 'maximum-extent' to override layer extents)

And if I define an envelope, I can only manage to have a cropped view of the earth with

mapnik.Envelope(-4500000, -4500000, 4500000, 4500000)

or a blank image with

mapnik.Envelope(-4600000, -4600000, 4600000, 4600000)

EDIT : I also tried specifying maximum-extent in the xml file but the result was the same as specifying the envelope in the python script.

What am I doing wrong ? Any hint will be appreciated :)

My configuration is :

WGS84 shapefile coastlines-split-4326.zip (http://openstreetmapdata.com/data/coastlines)
mapnik 2.2.0-pre
archlinux

This is the python script that renders the image :

#!/usr/bin/python2

import mapnik

image = 'world-xml.png'
map_config = 'world.xml'
m = mapnik.Map(1000, 1000)
mapnik.load_map(m, map_config)
bbox = mapnik.Envelope(-4500000, -4500000, 4500000, 4500000)
m.zoom_to_box(bbox)
mapnik.render_to_file(m, image)

And this is the xml config file :

<Map background-color="#ffffff" srs="+proj=ortho +lat_0=0 +lon_0=0 +ellps=WGS84 +units=m +x_0=0 +y_0=0 +no_defs">
        <Style name="projet-carte">
                <Rule>
                        <PolygonSymbolizer fill="#ffffff" />
                        <LineSymbolizer stroke="#000000" stroke-width="0.1" />
                </Rule>
        </Style>
        <Layer name="world" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
                <StyleName>projet-carte</StyleName>
                <Datasource>
                        <Parameter name="type">shape</Parameter>
                        <Parameter name="file">lines.shp</Parameter>
                </Datasource>
        </Layer>
</Map>


texts: Intermediate files disappear in ArcGIS 10.1

texts:
I am working with ArcGIS 10.1 creating some slope and aspect files. I do this:

DEM_park = ExtractByMask(mosaic_name, park_area)

and it creates the DEM_park file in a default gdb with a default name. I then export this file with:

output = 'dem_'+'{0}'.format(parkID)+'.tif'
arcpy.CopyRaster_management(DEM_park, output, "#","#","#","#")

which works fine. Then I want to keep working with the DEM_park file:

outSlope = Slope("DEM_park", "DEGREE")

But it says that the file does not exist. And indeed when I go the the default gdb I cannot find the file.

I can work with the .tif file I created, but still find it irritating that my other file just disappeared. Is there a way of avoiding this?

Additional Note: this happens on a script. If I execute the lines one by one in the Python window it can find the file and calculates the slope without problem.


texts:
Am trying to mosaic 11 scenes in ArcMap but my results are not smooth.

How can I mosaic and at the same time feather in ArcMap 10.1 to make sure the brightness is uniform.

I have EVI geotiff scenes


texts: Open source web map client with time sliders using WMS source

texts:
I have a web map service (WMS on geoserver). It has a current visualization which is decent. The data is divided into months currently using tabs. I would like to visualize the data using an open source based web map client with a time slider.

I have seen some examples around, but I am looking for a generic solution where I basically can change background layers (also WMS), source wms and relevant attributes. This to be able to experiment further with the solution. The wms time layer (wms-t) should not change the background wms-layers. They should remain controlled by a config-file or interactive selection interface not related to the time slider.

The current services can be accessed through the wms time parameter:

time=2014-1
time=2014-5
etc

It should be possible for the user to set min/max values and intervals, or the client should be able to guess this by itself based on the getcapabilities-info.

ESRI or other commercial (closed source) solutions are not possible options for this question.


texts: arcgis 10.2 javascript api proxy printtask error

texts:
when i invoke a printtask without the proxy, everything runs fine. But, when i use it, the invocation fails.

It seems that the cause is that the proxy is appending a double quote to the end of the url, so instead of sending

...Layout_Template=test

it sends

...Layout_Template=test"

Details:

Java proxy v 1.0


texts: Find or create XYZ tiles with only country names

texts:
I need XYZ tiles with only country names to use in a Leaflet website.
Do publicly available XYZ tiles with only country names exist? If so, where can one find such resources?

If not, is there such a way to generate XYZ tiles only with country names, controlling sizes and etc for the zoom level?

I would like to have something similar to Stamen Terrain Labels, but they contain city, continents and other information that are not country names.

Stamen Tiles are open source via Github, is it possible to get it into QGIS 3 or a similar visual software, in order to modify it to remove city names and just keep countries? If so, how?


texts: Limiting map tips popups by scale in QGIS?

texts: I have a points later and have enabled map tips to appear when I hover, as normal. However, I would like to prevent map tips balloons from appearing when zoomed out. The way I have found to accomplish this so far is to duplicate the points layer and set the point symbols to transparent, enable map tips in that layer, and disable in the original layer, and then limit the visibility of this new layer entirely. However, I'd like to think there is a more elegant way of limiting the popups from appearing when zoomed out.


texts: vb.net arcobjects check if topology has error

texts:
I follow the esri example for validating topology, is it possible to tell whether topology error exists? The code I use:

    ' Step 1 - access a topology
    Dim topoUiD As UID = New UIDClass
    topoUiD.Value = "esriEditorExt.TopologyExtension"
    Dim m_application As IApplication = My.ArcMap.Application
    Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD), ITopologyExtension)  'original code from esri
    Dim topology As ITopology

    topology = CType(topologyExt.CurrentTopology, ITopology)
    MsgBox("successfully access current topology")

    ' Step 2 - Finding topology errors
    Dim geoDS As IGeoDataset = CType(topology, IGeoDataset)
    Dim errorContainer As IErrorFeatureContainer = CType(topology, IErrorFeatureContainer)
    Dim eErrorFeat As IEnumTopologyErrorFeature
    eErrorFeat = errorContainer.ErrorFeaturesByRuleType(geoDS.SpatialReference, esriTopologyRuleType.esriTRTLineCoveredByLineClass, Nothing, True, False)
    Dim topoError As ITopologyErrorFeature

    ' check the eErrorFeat length to determine whether has error
    topoError = eErrorFeat.Next
    MsgBox("check topology errors.")

I am using vs 2012 and arcmap 10.2.


texts: I installed the plugin qgsaffine but cannot find the entry in the menu bar. Where do I find this. Can anyone here help me with this?


texts: How best to map lines based on multiple attributes?

texts:
I need to plot lines onto a map to display road quality. I'm already using colours to distinguish between different road conditions, but I also need to display another factor - the survey origin.

There are only two different origins to display, but I'm not sure what's the best way to do this. Dashed lines with gaps are not appropriate as surveys can legitimately contain gaps.

My current approach is to display the two types as follows:

Solid colour lines for one of the survey types.
For the second type, use a dashed line, but with the gaps drawn in a lighter shade of the same colour (shown below).

However it still seems a little confusing to see where surveys start and end.

Are there any alternatives or best practices to follow?


texts: How do I make a grid around all my points, and then convert it to a polygon

texts:
I want to make a grid around all my points, convert to a polygon and use this to clip the extent of my data. I am using Arc 10.2.

How do I make a grid and then convert it to a polygon?
Or is there a way to miss some of these steps?


texts:
I am serving several raster images with Geoserver using ImageMosaic and PostgreSQL to index data. The files are stored in the file system in GeoTiff format with EPSG:4326 and a pixel size of 2.7e-5 x 2.7e-5 arcsecond/pixel.

When I open this file in QGis with its vector shape definition, they match up perfectly in any OTF configuration, as it can be seen here:

The problem is that when the image gets served in Geoserver it seems to get skewed or shortened in the vertical axis, as seen here:

Here's another example. The raster image has a 1.5 pixel buffer that is why it should overflow the boundary, but the reprojection is just wrong.


texts:
I'm having difficulties to get the bounding box of a polygon. My code is the following:

def getBoundingBox(self):
    boundingBox = []
    polygon = self.getVectorLayer()
    vertices = polygon.geometry().asPolygon()

    xmax = vertices[0].x()
    xmin = vertices[0].x()
    ymax = vertices[0].y()
    ymin = vertices[0].y()

    for vertice in vertices:
        if vertice.x() > xmax:
            xmax = vertice.x()
        if vertice.x() < xmin:
            xmin = vertice.x()
        if vertice.y() > ymax:
            ymax = vertice.y()
        if vertice.y() < ymin:
            ymin = vertice.y()

    boundingBox.append(xmax)
    boundingBox.append(xmin)
    boundingBox.append(ymax)
    boundingBox.append(ymin)

    return boundingBox

The error message I get is

  AttributeError: 'QgsVectorLayer' object has no attribute 'geometry'

The function self.getVectorLayer() works perfectly fine for other tools I wrote, thus this shouldn't be the problem. With the line xmax = vertices[0].x()
I wanted to get the first vertice in the polygon, which probalby doesn't work that way, but it seems the problem starts even earlier.

Any ideas?


texts: Syntax error using layer.getFeatures PyQGIS

texts:
I am simply trying to run a command in QGIS to extract the coordinates of a point layer, but I get a syntax error when trying to do so (see code below). I can not see where the issue is - any suggestions?

for f in layer.getFeatures():
  geom = f.geometry()
  print geom.asPoint()

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 63, in runsource
    code = self.compile(source, filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 168, in __call__
    return _maybe_compile(self.compiler, source, filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 99, in _maybe_compile
    raise err1
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 87, in _maybe_compile
    code1 = compiler(source + "\n", filename, symbol)
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 133, in __call__
    codeob = compile(source, filename, symbol, self.flags, 1)
  File "<input>", line 3
    print geom.asPoint()
             ^
SyntaxError: invalid syntax


texts: Limiting decimal places of coordinate values in GEOJSON file using QGIS?

texts:
Using "QGIS 2.0.1-Dufour", I try to convert a vector layer into GEOJSON file using "Save As" sub menu. I get the GEOJSON  file, but the coordinates of points have values up to 15 decimal places.

How can I round off the values at the desired number of decimal places?


texts: Where to find surface level, monthly wind speed and direction data

texts: I'm looking for monthly averages for wind speed and wind direction at surface level for South America, preferably (but not exclusively) as GIS layers.


texts: Dealing With ArcPad and Geodatabase Domains

texts: I created a polygon feature class through a file geodatabase. In the geodatabase i created domains and gave them drop downs. When i go to create a feature class in editor in ArcGIS, in the attribute i am allowed to see dropdowns. However when i use ArcPad i get a empty text box rather than a drop down. What am I doing wrong?


texts: How to avoid hangup when clipping a large polygon in QGIS

texts:
I am working with QGIS 2.12.3 Lyon. I have a shapefile which contains slope data for a large territory (layer - PenteSZ32). I want to lighten the weight of my shapefile by clipping it to a narrower surface that correspond to my study area. My study area correspond to a ripperian buffer of 1000 meters (layer - Donut_buffer_1000m).

In order to clip my polygon, i used the Vector > Geoprocessing Tools > Clipper as shown in the following image. This request worked previously with other layers. However, it doensn't work in this case. It bugs at 12 %. Both layer are in the same SCR.

Can you help me figure out what is the problem with my request ? Alternatively, do you have any other way to proceed you could suggest me ?


texts: Area and distance calculation/measurement unexpected after projecting?

texts:
I have converted a GIS shapefile of Botswana from the geographic coordinate system GCS_WGS_1984 to the projected coordinate system Africa_Albers_Equal_Area_Conic to be able to calculate the area of the polygons of the shapefile (and later to calculate buffers etc.). But the scale of the projected shapefile is completely wrong, which skews the area (and distance) calculation. For example, when I use the ruler to calculate the distance between two points, a ~113km distance in reality is measured as a 0,001087 km distance on the data frame, thus about 100 000 times too small...

But when I do the same distance measurement with the ruler when the shapefile is under the geographic coordinate system GCS_WGS_1984, there is no problem.

Does someone have an idea how to fix this problem, and end up with a correct scale when I'm using the projected coordinate system?


texts:
I am using slope function from geoHMS Toolkit in arcMap 10.3. After the execution of the function the river attribute table is showing 0 in all the fields of slope,ElevUP  and Elev down. I am projecting the river with same coordinates as rawDem.

I am attaching the attribute table and projection of the river(stream).


texts: Discrepancy between DXF and Shapefile Polygons

texts:
I have created a shapefile based on a dxf.

the areas are as follows:
DXF: 51216417 square metres = 5121.64 Hectares = 51.21 square Kilometres.
GIS: 50.998 square kilometres.

What could have caused this?

Running Autocad 2013 and QGIS 1.8


texts:
I have a raster DEM and a polygon covering a part of it. I need an area-elevation distribution (i.e. how many m2 are in which elevation interval). 20m intervals will be fine.
I can't figure it out.


texts: Changing symbology/color in kml when exported from ArcGIS/QGIS?

texts:
I have marked some points and used symbology in ArcGIS/QGIS but when I export the file to kml, the color/symbology changes and only single color/symbol remains.

Is there any possible way that the symbol/color that I use in ArcGIS or QGIS may retain in kml as well?


texts: What is the 'Layer' box in 'Save Vector Layer' good for?

texts: In QGIS, in the 'save vector layer as' window, under OGR creation options, what is the 'Layer' box for?


texts:
First of all I'm a complete noob when it comes to QGIS so maybe thats the reason why I cannot find anything on this on google (or is this due to the 2.8.1 compatebility issue?).

My issue is that when I create the heatmap the location of the heatmap is off compared to the data. I've checked the CRS and all match (project CRS, data CRS AND the heatmap CRS). I've redone several tutorials several times and started several times from scratch but somehow it keeps giving me the data on the wrong location :-(.


texts: Generating CityGML file with number of storey attribute values from shp file?

texts:
I am trying to convert shp file to cityGML. In shapefile I have building footprints with attributes "storey" as number of storeys. I want to match this attribute to number of storeys attributes of cityGML buildings class.

How can I make this?


texts: Querying Raster Dataset by Location to Study Area

texts:
I have downloaded a very handy raster dataset generated by the Nature Conservancy for identifying terrestrial habitat cover types in the north east.  Link to the project info and a link containing a download of the raster dataset (right hand side of page) is here:

https://www.conservationgateway.org/ConservationByGeography/NorthAmerica/UnitedStates/edc/reportsdata/terrestrial/habitatmap/Pages/default.aspx

I'm not terribly experienced when it comes to working with raster files. All I intended to do with this data was query which habitat types occur within a quarter mile of a series of reservoirs in New York. However, select by location is not an option with rasters. I've tried converting the raster to polygon but not only are the polygon boundaries slightly off, there is no information retained in the resulting attribute table that would let me join the original dataset info to the resulting shapefile. Does anybody know how I can go about inventorying the habitats in my study area?

I have attempted to research a solution however I've come up short. I have access to Spatial Analyst and an Advanced license.


texts: Seeking official topographic maps of Cuba?

texts: Any ideas as to where I might find official topo maps of Cuba at scales of !:50K or larger?  NGIA won't release the ones they have, and similarly they can't be purchased in Cuba.


texts:
I would like to use a basemap from the Ancient World Mapping Centre, available as maptiles via Mapbox, as basemap in my ArcMap desktop version 10.2.2.

I tried to import it as a new WMTS server with the following url:

http://api.tiles.mapbox.com/v4/isawnyu.map-knmctlkh/{z}/{x}/{y}.png?access_token=YOUR_TOKEN_HERE

I replaced the 'YOUR_TOKEN_HERE' by my token from Mapbox, but ArcMap did not want to create the WMTS.

I tried it in QGIS using the Tile Layer Plugin and everything worked perfect.

I prefer working with ArcMap, so I would like to know how I can import these tiles from Mapbox in my ArcMap desktop.


texts: Get a window from a raster in rasterio using coordinates instead of row/column offsets, width and height

texts:
I am trying to open a windowed dataset in rasterio, but I have a bounding box to work with. I thought it would be easy to do this but it seems that it involves some complicated process with a "Mixin", something I am having a very hard time understanding, even after reading the thread here.

The documentation says:

  A subclass with this mixin MUST provide the following properties:
  transform, height and width

But I don't know how I can implement these properties if I don't know ahead of time what the height and width are going to be, and I'm not even sure that I am implementing the code properly anyway. This is what I have so far:

import rasterio
from rasterio.windows import WindowMethodsMixin, Window
from rasterio.enums import Resampling

class MyWindow(WindowMethodsMixin, Window):
    pass

with rasterio.open("flask/Docker/app/dem/slope_sm.tif") as src:
    rst1 = src.read(1, window=MyWindow.window(...))

I know this is hardly complete but I am really confused about how to proceed.

My IDE tells me that the call to MyWindow.window() takes the following parameters:

self
left
bottom
right
top
precision (optional)

But I don't know what to pass to the "self" parameter. The rest of the parameters are understandable - they are the edges of my bounding box from which I am constructing the window.

The documentation doesn't supply any working example of this, which I find odd because this seems to me like a pretty typical operation.

Can somebody please show me how to do this, or at least explain what I am missing here?


texts: Getting trouble with prefix clause for an ArcObjects Java query filter

texts:
I'm using ArcObjects/Java for ArcGIS 10.2.2 and I'm trying to get unique values of a field for a feature class.

Here's my current code:

List<String> values = new ArrayList<String>();

int fieldId = featureClass.findField(fieldName);

IQueryFilter queryFilter = new QueryFilter();
queryFilter.setWhereClause(whereClause);

((IQueryFilterDefinition2)queryFilter).setPrefixClause("DISTINCT " + fieldName);

IFeatureCursor cursor = featureClass.search(queryFilter, false);
IFeature feature = null;

while((feature = cursor.nextFeature()) != null) {
    values.add((String)feature.getValue(fieldId));
}

Cleaner.release(cursor);
return values;

In theory, by using the 'prefix clause', the SQL statement should be SELECT DISTINCT field FROM ... WHERE ... but when I execute this code, I'm getting an:

  AutomationException : An invalid SQL statement was used. [SELECT
  DISTINCT CLASS_ADM * FROM TRONCON_ROUTE]

.

The problem is this * after the DISTINCT clause that shouldn't be here.

I tried to set subFields to nothing like:

queryFilter.setSubFields("");

or

queryFilter.setSubFields(null);

but this doesn't solve the problem.

I also tried to use IDataStatistics.getUniqueValues but it returns a deprecated type IEnumVARIANT that I don't know how to use..

Any idea to get that working?


texts:
I have a geotiff that is around 26MB on my hard drive. I would like to use leaflet to plot it within a shiny application.
However, it just doesn't load. When I use mapview just to plot it locally it also takes some seconds, but eventually it will be plotted. To my knowledge it's better to use leaflet in a shiny-app. So I don't know exactly what to do now.
Here is some code, that produces a raster of more or less similar size
r = raster()
res(r) = 0.1
r[] = runif(ncell(r))

# plot it with mapview
mapview::mapview(r)

# plot it with leaflet
leaflet() %>% addTiles() %>% addRasterImage(r)


texts: How can I automatically make new computings after I have modified some value in a field in ArcGIS Desktop?

texts:
I'm a computer programmer and I have no idea about geographic things.
I have got a long algorithm with a lot of computings, computed from a previous computed field. There are 7-8 computed fields and they refer to each other, so it's a bit complex algorithm. There is a chance, that the original values in the beginning of the long algorithms will change (some empty columns will have values and the algorithm needs to compute from them). How can I make this to be as automatic as it can be? I think there is other solution than run the python scripts on the fields again and again.


texts: Showing information of several layers of GeoJSON request to GeoServer?

texts:
Always enters the error function and shows me the alert 'what happen?'.

This is the piece of code that causes problems:

map.addEventListener('click',onMapClick);
function onMapClick(e)
{
    var latlngStr = '('+ e.latlng.lat.toFixed(3) + ", " + e.latlng.lng.toFixed(3)+')';
    var BBOX = map.getBounds()._southWest.lng + ", " + map.getBounds()._southWest.lat + ", " + map.getBounds()._northEast.lng + ", " + map.getBounds()._northEast.lat;
    var WIDTH = map.getSize().x;
    var HEIGHT = map.getSize().y;
    var X = map.layerPointToContainerPoint(e.layerPoint).x.toFixed(0);
    var Y = map.layerPointToContainerPoint(e.layerPoint).y.toFixed(0);
    URL = 'http://192.168.60.55:9090/geoserver/ems/wms?request=GetFeatureInfo&service=WMS&version=1.1.1&info_format=application/json&exceptions=application/vnd.ogc.se_xml&layers=ems:tblsevilla,ems:tblmanz,ems:tblcarre&styles=&srs=EPSG%3A4326&format=image%2Fpng&bbox='+ BBOX + '&width='+ WIDTH +'&height='+ HEIGHT + '&query_layers=ems:tblsevilla,ems:tblmanz,ems:tblcarre&feature_count=20&x='+ X +'&y='+ Y;

    function onEachFeature(feature, layer) {
        if (feature.properties && feature.properties.PERIMETRO)
        layer.bindPopup("Perimetro: "+feature.properties.PERIMETRO);
        }
    }

    $.ajax({
        type: 'POST',
        url: URL,
        dataType : 'json',
        success: function (response) {
            alert("2");
            L.geoJson(response, {
                onEachFeature: onEachFeature
            }).addTo(map);
        },
        error: function (req, status, error) {
            alert('what happen?');
        }
    });
}

I'm a Leaflet and GeoServer beginner.


texts:
Is it possible to generate some ol3 object in node server then put them into response?
I want to do that for line string


texts: Coordinate system and projection for georeferencing Google Maps image in ArcGIS Desktop?

texts:
I have saved Google Maps image (not Google Earth) of area at latlon(17.692437,83.218689) by some method in jpeg format. It is an exact snapshot of the Google Maps. Now I need to geo-reference it to use it as base map for my digitization. What projection or coordinate system shall I use? Some says its WGS-84 Mercator projection. But I see no Mercator projection in GCS list in ArcGIS. Please help! I am confused.
I want to georeference it by assigning exact latlon values to the points as in google maps.
Extend : Now I got another question. Is Google Maps in Projected Coordinate system or Geographic Coordinate system. Though it displays LatLon values, I think it is projected coordinate system. Am I right?
Extend
I set the projection of ArcMap Dataframe to Projected Coordinate System ( Web Mercator) and changed the display values from Meters to Decimal Degrees. Then I added the Google Map with no geo-referenced data to dataframe. Then I took some common control points in my ArcMap image and Google Maps, set the Latitude and Longitude values from Google Maps in ArcMap image points, After assigning latlon values to 4 points, I click on update georeferencing. Then my image is transformed very strange. Is the procedure I did correct?


texts: How to use QGIS ComboManager inside my plugin?

texts:
Could someone provide some detailed instructions (for very new beginners!!) on how to integrate combomanager in another plugin?

What I have tried (maybe wrong way of doing it!!...):
First I downloaded Qgis Combo Manager from http://3nids.github.io/qgiscombomanager/ and extracted the zip folder.

Then I copied
composercombo.py, expressionfieldcombo.py, fieldcombo.py, layercombo.py and optiondictionary.py in my class directory.
Then I included in the import section of my module' s python file
the following:

from qgis.core import *
from layercombo import VectorLayerCombo, LayerCombo
from fieldcombo import FieldCombo
from expressionfieldcombo import ExpressionFieldCombo

Finally, under run(self)

line186    def run(self)
line187    self.dlg.comboBox = VectorLayerCombo(self.dlg.comboBox)
line188    self.dlg.comboBox_2 = FieldCombo(self.dlg.comboBox_2, self.dlg.comboBox)

So when I run my plugin for the first time it works. Then, if I run it again just after the first time I get this error message:

Traceback (most recent call last):
  File "C:/Users/Eleftheria/.qgis2/python/plugins\MyClass\my_module.py",line 188, in run
self.dlg.comboBox = VectorLayerCombo(self.dlg.comboBox)<br/>
  File "C:/Users/Eleftheria/.qgis2/python/plugins\MyClass\layercombo.py",line 187, in __init__
 LayerCombo.__init__(self, widget, initLayer, options, QgsMapLayer.VectorLayer)
  File "C:/Users/Eleftheria/.qgis2/python/plugins\MyClass\layercombo.py",line 62, in __init__
    self.widget.currentIndexChanged.connect(self.__currentIndexChanged)
AttributeError: 'VectorLayerCombo' object has no attribute 'currentIndexChanged'
Python Έκδοση:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
Έκδοση Qgis:
2.10.1-Pisa Pisa, d20c5b7

In line 62 there is this:

line62    if qVersion() > '4.3.3':
line63        QCoreApplication.installTranslator(self.translator)


texts: How to add non spatial sqlite table to map using python in qgis?

texts:
Is it possible to load pure sqlite 2.0 (without spatial extension) table into qgis using python?  I would like to get the same result as in the case of using Add Vector Layer -> All files.
I have already modified this example http://www.qgis.org/pyqgis-cookbook/loadlayer.html
but it doesn't work in my case.
Can someone help me?
I use in  the python console this code:

>>> db = r'c:\Program Files (x86)\Quantum GIS Lisboa\data\base.db'
>>> uri = QgsDataSourceURI()
>>> uri.setDatabase(db)
>>> schema = ''
>>> table = 'table'
>>> geom_column = ''
>>> uri.setDataSource(schema, table, geom_column)
>>> display_name = 'table'
>>> vlayer = QgsVectorLayer(uri.uri(), display_name, 'sqlite')
>>> QgsMapLayerRegistry.instance().addMapLayer(vlayer)
<qgis.core.QgsVectorLayer object at 0x0CFA08A0>

i also tried with 'spatialite' in vlayer. When I use Add vector layer i can see result table in the TOC. Using this code I can't display table in table of contents.  I think the reason of my trouble is trivial, but I have no idea what I should do to see my table


texts: Application Templates for ArcGIS Server?

texts:
Esri provide a lot of excellent template applications for us in conjunction with web maps.

These templates are specifically designed for using with your resources on ArcGIS Online as opposed to ArcGIS Server.  (i.e. you supply a web map ID, as opposed to URLs to your GIS Services Directory).

In particular I found this awesome template, that ticks all the boxes for me:

Mobile support (e.g media queries, responsive design)
Share to social media
Nice looking layout
Built in Geocoder
etc. etc.

I have been looking at the code behind it, and was hoping I could reverse engineer it to work with my ArcGIS Server REST endpoints instead.  However, does not look like an easy task (im new to Javascript) and believe I would spend more time trying to do this Vs building something from scratch.

Can anyone recommend some template applications to use as a starting point for developing an ArcGIS Javascript API app to use against ArcGIS for Server?

I have already put together the bones of my application in a basic dojo template, but its the stuff mentioned above that I need to also include in my web application.

This will then allow me to focus my development efforts on the actual Esri Javascript API logic, without needing to spend a lot of time having to do a 101 on other aspects of web development.

I understand there are a number of web frameworks out there (e.g Bootstrap) that might assist in some of this, but I was kind of hoping for something a bit more established and aimed at Esri Javascript API.  However, would still be interested in answers that do leverage other web framework templates as a starting point.


texts: Is the QGIS Development nightly repo nightly?

texts:
I've recently installed QGIS 3 (2.99) on Ubuntu 16.04 Xenial using the Nightly build for Debian and Ubuntu  following http://www.qgis.org/en/site/forusers/alldownloads.html . The installation was successful, and I was able to use the Development version of QGIS3 for most of the functions so far.

However, I've noticed that my version does not seem to be up-to-date with the nightly source code??, despite the fact that I've updated my Ubuntu several times.

I noticed this because I encountered the schema "g" does not exist bug https://issues.qgis.org/issues/15910, that was reported fixed a couple of months ago. Check my QGIS code revision, it seems to be something committed half a year ago

Is the Nightly build for Debian and Ubuntu really nightly w.r.t QGIS code revisions?

Also, is there a way to ensure that my repo is updated? I'm thinking about removing QGIS entirely, and then reinstall. But I'd better ask here first.


texts: ArcPy script runs longer processing time in scheduled task than in ModelBuilder?

texts:
I have a Python script that runs ArcPy and a model that is imported with a Toolbox.
When directly running the model on ArcMap, it takes 4 plus minutes. It takes about the same time if running in the Windows Command Prompt window. E.g.
Python.exe "{SCRIPT_PATH}"

However, if I run the script in a Windows Scheduled Task, it takes 8 minutes to complete.
One of the processes in the model that makes a big difference is "Tabulate Area". It takes 14 seconds in ArcMap while 44 seconds when run as Scheduled Task!
Running in command prompt and running the scheduled task are done on the same machine.
My environment: ArcMap 10.1, Windows Server 2008 R2 Enterprise, Python 2.7.2.
Why is there this performance difference?


texts: Set projection from user input with Arcpy

texts:
This seems like to should be an easy thing to do, but I can't seem to make it work.

arcpy.env.workspace = arcpy.GetParameterAsText(0)
output_working = arcpy.GetParameterAsText(1)
PLSS_Intersected = arcpy.GetParameterAsText(2)
spatial_ref = arcpy.GetParameterAsText(3)

outCS = arcpy.SpatialReference("spatial_ref")
outfc = os.path.join(output_working, "PLSS_Intersected")
arcpy.Project_management(PLSS_Intersected, outfc, outCS)

EDIT:

How is the error massage that I receive. The very top is what prints out when I printed outCS.


texts: Using Python function for cost calculation in GIS

texts:
I want to calculate the cost of two raster layers (ramdom, simulated) composed of values 1 and two, the first created randomly and the second by optimization using two associated cost layers. The code looks somewhat like below. However, I cannot make the function to read my second raster.

driver1 = gdal.GetDriverByName ('RST')
driver1.Register ()
ready = gdal.Open ('ramdom.rst', gdal.GDT_Byte)
cols = ready.RasterXSize
rows = ready.RasterYSize
bands = ready.RasterCount
geotransform = ready.GetGeoTransform()
read_ramdom = ready.ReadAsArray ()

global sum1, sum2, value1, value2
sum1, sum2, value1, value2 = 0, 0, 0, 0
classes = 2
global ssum
def cost(something):

          global ssum
          global sum1, sum2, value1, value2
          global classes
          while classes >=0:
             if something == something:
              something = ready.ReadAsArray()
             elif something == something2:
              something = simul_read.ReadAsArray()

              for i in range (0, rows):
                   for j in range (0, cols):
                      numb =  something [i, j]
                      if numb == 1:
                       one = names [numb -1]
                       ds1 = gdal.Open (one, gdal.GDT_Float32)
                       data1 = ds1.ReadAsArray (j, i, 1, 1)
                       value1 = data1 [0, 0]
                       sum1 = sum1+value1

                      if numb == 2:
                       two = names [numb -1]
                       ds2 = gdal.Open (two, gdal.GDT_Float32)
                       data2 = ds2.ReadAsArray (j, i, 1, 1)
                       value2 = data2 [0, 0]
                       sum2 = sum2+value2
                       classes = classes - 1

          ssum = sum1 + sum2
          return ssum

simul_read = gdal.Open ('simulated.rst', gdal.GDT_Byte)
read_simulated = simul_read.ReadAsArray ()

A = cost(something)
B = cost(something2)
C = cost(something) - cost(something2)

print ' The initial COST was :','   ', A, '\n'
print ' The final COST is :','   ', B, '\n'
print ' The gain is :', C , '\n'


texts: Setting the Centroid of a Set of Polygons

texts:
I'm trying to shift a set of rectangle polygons over an area of interest (AOI).  I thought it would be as easy as setting the SHAPE@TRUECENTROID of the polygon area of interest to the shapefile/featureclass of the rectangles. When multiple polygons are selected in Arcmap, the centroid of the collection of polygons appears.

When attempting to pass the shapefile in to an update cursor, each polygon overlaps over the AOI coordinates instead.  Is there a simple, Arcpy, way of sifting the entire set of polygons in a shapefile to be centered over said area of interest?  Or put another way, maintain the spatial relationship each rectangle has with one another while shifting the entire set over the AOI.

I should add that I am first rotating the polygons with Arcpy.  When running the script it doesn't keep the rotation expected.  I would expect setting the centroid would merely shift, not rotate.

EDIT.
Here is the rotate code (this was adapted from another post on here dealing with rotating polygons)

def RotateAxis(AnchorX,AnchorY,inputx,inputy, rAngle):
        x = inputx - AnchorX
        y = inputy - AnchorY
        resultx = ((x * math.cos(rAngle)) + (y * math.sin(rAngle))) + AnchorX
        resulty = ((x * math.sin(rAngle)) - (y * math.cos(rAngle))) + AnchorY
        return (resultx,resulty)

I update the geometry from this then finally after each polygon has been rotated I attempt to perform the shift.

XC and YC are the anchor points which are the SHAPE@XY from my AOI.

Updated Code to reflect bottom picture.

memFC = "in_memory\\dissolvedPoly"
arcpy.Dissolve_management(inFC, memFC,"PageNumber")

with arcpy.da.SearchCursor(memFC, ["SHAPE@XY"]) as oFeature:
        for row in oFeature:
            memX, memY = row[0]

xShift = float(XC) - float(memX)
yShift = float(YC) - float(memY)
print "Shift x " + str(xShift) + " ; Shift y " + str(yShift)

with arcpy.da.UpdateCursor(inFC, ["SHAPE@XY"]) as cfeature:
        for row in cfeature:
            print "Shifting polygons..."
            print ""
            print "x = " + str(x), " y = " + str(y)

            cfeature.updateRow([[row[0][0] + (xShift), row[0][1] + (yShift)]])

Before Shifting Polygons(What I start with)

After Shifting Polygons (What I Want)

What I actually get.

EDIT 2:

What I get with the shift x, shift y code.


texts:
I have an Image with the following metadata:

Projection Information :

    Projection Zone  : 0
    Spheroid Name   : Modified Everest
    Georeferenced to : RSO

Map Information :

    Upper Left center X : 598398.000000
    Upper Left center Y : 506921.000000
    Lower Right center X: 622798.000000
    Lower Right center Y: 473521.000000
    Pixel X size : 20.000000
    Pixel Y size : 20.000000

The Geographic Coordinate System the image uses is:

Geographic Coordinate System: GCS_Kertau_RSO
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_Kertau_RSO
 Spheroid: Everest_1830_RSO_1969
  Semimajor Axis: 6377295.664
  Semiminor Axis: 6356094.667915204
  Inverse Flattening: 300.8017

I need to project this image in the WGS_1984_UTM_Zone_48N coordinate system. I am using the ArcGIS-Desktop package. To achieve this, I am using the Project Raster tool under the Data Management Tools to convert the projection. The problem is that even after the projection takes place, the projected image does not get properly aligned to its place and stays at the same spot (over the original image).

On checking the layer properties of the projected image, it shows the layer is spatially referenced to the WGS_1984_UTM_Zone_48N system, but unfortunately, it is not displayed at its correct position (which can be checked by adding a BaseMap from ArcGIS online resources).

Need suggestions/ alternative solutions to the problem.


texts:
Is it somehow possible to change the buffer style of the build in GeoTools buffer() function? Or is there a User Definded Function for that? I just saw a different buffer function where one could define the number of segments.

I need something similar to the options QGIS has. I want to define the "End cap style" and the "join style". At the moment GeoTools creates the buffers with round corners. But when I for example buffer a triangle I want it to remain a triangle just with a greater area.


texts: I have loaded a shapefile into QGIS and I now want to export it into our Oracle Database and add SDO_GEOMETRY. Is it possible to do this directly in QGIS? If not, what opensource options do I have to perform this task?


texts: Running basic JavaScript in KML to be viewed through OpenLayers API

texts:
I am trying to build an application in which points can be loaded on OSM map using KML called through OpenLayers with a balloon popup on each point. I am using the sample Sundials application: http://jorix.github.io/OL-FeaturePopups/examples/sundials-FP.html

I want to put a very simple album in the balloon using which the user can see more than one images within the balloon itself. I am sure this cannot be achieved without running some JS within the application. I thought I had a chance when I found the KMLs at: http://kml-samples.googlecode.com/svn/trunk/kml/javascript/

The scripts in these KMLs, when opened in Google Earth, works just fine but don't when called through OpenLayers. Upon digging in the developer tools it was revealed that the script part of the CDATA was not accessed and hence the function was declared undefined.

What would be a workaround for this? Any other way to view an album in a map popup?


texts: Convert coordinate reference system of GeoJSON using ogr2ogr with same file name?

texts:
I am currently looking to use ogr2ogr to convert the coordinate reference system of a GeoJSON file. I want the output file to have the same name as the input file. This is on a Linux PC running Ubuntu 20.04 LTS on Oracle Cloud.
This is the code I tried first with the output below
ogr2ogr -s_srs EPSG:32614 -t_srs EPSG:4326 -f GEOJSON plot_boundary.geojson plot_boundary.geojson

ERROR 1: Source and destination datasets must be different in non-update mode.

This is the code I tried second with the output below
ogr2ogr -update -s_srs EPSG:32614 -t_srs EPSG:4326 -f GEOJSON plot_boundary.geojson plot_boundary.geojson

ERROR 1: Layer plot_boundary already exists, and -append not specified.
        Consider using -append, or -overwrite.
ERROR 1: Terminating translation prematurely after failed
translation of layer plot_boundary (use -skipfailures to skip errors)

This is the code I tried third with the output below
ogr2ogr -overwrite -s_srs EPSG:32614 -t_srs EPSG:4326 -f GEOJSON plot_boundary.geojson plot_boundary.geojson

ERROR 5: -nln name must be specified combined with a single source layer name,
or a -sql statement, and name must be different from an existing layer.


texts:
I will save my project into Map Package so I can share with another user. Using ArcGIS I can easily make Map Package as .mxd but I don't know how to do that on QGIS.

So how to make Map Package on QGIS?


texts: Visualizing updating lines of bearing from point?

texts:
I'm working on a project where I need to visualize a constantly updating line of bearing from a point, based on degree of rotation around that point. My plan is to design a line feature with an attribute table that links to a CSV which is being fed this rotation data every few seconds. I'll run a python script that creates a new line feature each time the CSV is updated, visualizing a max of 10 lines, removing the oldest line per update.

My big hurdle right now is figuring out how to make the generated lines rotate around the relevant point through maybe some sort of snapping feature.

Is this possible?

I feel like there's a much easier way to do this than the method I'm proposing, but I can't think of what it might be. This is ultimately going to be ported to a separate window to a UI outside of QGIS, so I'm worried I'm going about this the wrong way.

I've come across some similar questions that touch on what I'm getting at, but nothing that's truly addressed this specific need. I've also tried working with the Distance and Azimuth plugin, but it works in a way that wouldn't cooperate with my aforementioned method. I need the lines to snap to an independent point feature with separate attributes. Or maybe I'm wrong and this can be done.


texts:
I'm writing a plugin for QGIS in python

in the QGIS API doc, I can see the QgsDiagramRendererV2 class, under the qgs.core section
http://www.qgis.org/api/classQgsDiagramRendererV2.html

What I'm trying to do, is simply create diagram in a QgsVectorLayer.

so I try to do this :

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *

myLayer = QgsVectorLayer("some/url", 'label', "ogr")
myRenderer = QgsDiagramRendererV2()
myLayer.setDiagramRenderer(myRenderer)
QgsMapLayerRegistry.instance().addMapLayer(myLayer)

And then I get an error saying that the QgsDiagramRendererV2 class doesn't exist

But it exist in the doc

What can I do to instantiate a DiagramRender ?


texts:
I'm looking for some resources for learning PyQGIS.

It would be interesting having a collection of books or websites that provide some practical examples for learning the syntax or accomplishing specific tasks.

Ideally, these resources should give a general guidance for both beginner and experienced users.

Where to find QGIS tutorials and web resources? is a very similar question, but it gives help for learning QGIS, and not specifically PyQGIS (in fact, it doesn't have the PyQGIS tag).

Any help?


texts: Scatter plot of two columns for a FeatureCollection in Google Earth Engine?

texts:
I have some SNPP-FIRMS data in a FeatureCollection format. I would like to make a scatter chart of two of its columns (FRP Value, and Date). I tried but it did not give me a correct chart which does not show anything (see the screenshot).

)

I do not know why! I am not sure if it’s because of the format of the values of two columns (FRP: Float, and Date: Long).

Here is my code:

var viirs = ee.FeatureCollection('VIIRS_4Mar_21Oct_2018_T')

var Viirs_FRP = viirs.select('FRP');

var Viirs_Date = viirs.select('ACQ_DATE');

var chart_s = ui.Chart.feature.byProperty ({
     features: viirs,
     xProperties: ['Viirs_FRP'],
     seriesProperty:'Viirs_Date'
})
    .setChartType('ScatterChart')
    .setOptions({
      title: 'Fire Radiative Power (MW) per time (VIIRS)',
      hAxis: { title: 'DATE'},
      vAxis: { title: 'FRP VALUE CHANGE' },
      lineWidth: 1,
      pointSize: 3
});
print(chart_s);


texts:
I have a table on PostgreSQL 9.2 ("POSTGIS="2.0.1 r9979" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.8.0""). This is it:

CREATE TABLE nearest_point (
  id_nearest_point serial NOT NULL,
  id_iri_radartec integer,
  id_tramo integer,
  geom geometry(GeometryZ,4326),
  CONSTRAINT pk_nearest_point PRIMARY KEY (id_nearest_point)
);

I need a geometry that represents the point set union of the Geometries. Easy.

SELECT id_tramo
   ,ST_Multi(ST_Union(np.geom)) AS geom
  FROM nearest_point AS np
GROUP BY id_tramo;

But I get this:

"ST_Union GROUP BY" points in pink.

And the other points?

The arrow indicates where are the points of the table nearest_point.
Where is my error?


texts: Get already used shortcuts in QGIS plugin

texts:
I have a c++ plugin in QGIS which contains a treewidget. Now I want to get the keyPress-events. Until now I get this with the eventfilter described here

My problem is, I want to use shortcuts like Ctrl-N to create a new treewidgetitem. But this shortcut is used by QGIS self for creating a new document and it doesn't goes through my eventfilter. Knows someone a trick to get the event, if the focus is at my plugin-treewidget?


texts:
I have a shapefile with one layer. I can open it with QGIS. Now I want to convert this shapefile or layer to Spatialite database.

How can I do it? Is it possible to do with QGIS?


texts: Projected coordinates for Well-known ID 5232 cannot be found in ArcMap?

texts:
Projected coordinates for SLD99(WKID- 5232) cannot be found in ArcMap Coordinate system.
 Is there any built in procedure to find that original coordinate system in ArcMap 10.3? If not what is the procedure to define that coordinate system.


texts:
My attribute table (as in the attached screenshot) shows that there is only 1 element and the rest of the 26 items is filtered out.
How can I turn off this filter?


texts: Percentage of polygon at, or less than, specific height supplied by raster

texts:
I am using QGIS to process some LIDAR data. I need to be able to define the area of a pond with respect to its depth (i.e it is 100m^2 at 32mAOD and 500m^2 at 33mAOD).

Is there a way of finding the percentage (or area) of a polygon with its height at less than or equal to a specified value (the height being supplied from the LIDAR raster)?


texts: Why can't gdal save CRS information to GeoJSON?

texts:
gdal fails to add CRS definitions to GeoJSON files in some cases.

Try this, for example, from any
linux terminal that has gdal installed:

# download the data
wget "http://giscience.tw.rpi.edu/rpi_group_projects/watershedders/data/gas_station.geojson"

# convert the downloaded file to shapefile format
ogr2ogr -f "ESRI Shapefile" gas_station.shp gas_station.geojson

# convert the shapefile back into geojson
ogr2ogr -f "GeoJSON" gas_station2.geojson gas_station.shp

Both operations work fine, but the new geojson file
is missing a CRS. Here are the top 4 lines of
gas_station2.geojson:

{
"type": "FeatureCollection",

"features": [

Line 3, usually reserved for the CRS, is empty.
Yet the correct CRS is saved in the Shapefile's .proj file:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Note that gdal keeps CRS data when converting from GeoJSON to GeoJSON:

ogr2ogr -f "GeoJSON" gas_station3.geojson gas_station.geojson

# top lines of gas_station3:
{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

This "crs" member is the same as the original.

Is this a bug in gdal?

A similar issue was discussed
on the R-sig-geo email list.

A related bug report was opened 2 years ago, but this is now closed: http://trac.osgeo.org/gdal/ticket/4995

I've tried this on gdal 1.10.1 and 1.11.1.


texts: Gdal not working with Oracle Spatial's synonyms

texts:
I managed to get the OCI driver to work on both oracle table and views but I can't get it to work with synonyms.

Is this a limitation or am I doing something wrong?

Here's a sample of the command I'm using:

"C:\Program Files\GDAL\ogr2ogr.exe" -f "ESRI Shapefile" "test_shape.shp" "OCI:<schema>/<password>@(DESCRIPTION = (ADDRESS_LIST = (ADDRESS = (PROTOCOL = TCP)(HOST = localhost)(PORT = 1521)))(CONNECT_DATA = (SID =orcl))):MY_LOCAL_SYNONYM

I get the following error:

  ERROR 1: ORA-04043: object "MY_LOCAL_SYNONYM" does not exist

If I log into SQL Developer with the same schema/password and try to perform a select it works perfectly.

SELECT * FROM MY_LOCAL_SYNONYM;

SELECT * FROM <schema>.MY_LOCAL_SYNONYM;


texts:
Good morning,
I have a .csv file with more than 2000 rows that I would like to upload to qgis as point vector. When I add the text layer it loads only 1 value.
I have selected as X field the X coordinate and as Y field the Y coordinate.


texts:
i want to export a qgis legend with specified width , height and dpi. but it crops the content
layout = QgsPrintLayout( QgsProject.instance() )
layout.initializeDefaults()
layout.pageCollection().pages()[0].setPageSize(QgsLayoutSize( width, height, QgsUnitTypes.LayoutInches))
manager = QgsProject.instance().layoutManager()
manager.addLayout(layout)

legend = QgsLayoutItemLegend(layout)
legend.setTitle("title")

layerTree = QgsLayerTree()
for layer in layers:
    layerTree.addLayer(layer)

legend.model().setRootGroup(layerTree)
layout.addLayoutItem(legend)

exporter = QgsLayoutExporter(layout)
exporter_settings = QgsLayoutExporter.ImageExportSettings()
exporter_settings.dpi = 256
inch_to_pixel = 96
exporter_settings.imageSize = QSize( inch_to_pixel * width , inch_to_pixel * height )
exporter.exportToImage(path, exporter_settings)

when i remove layout.pageCollection().pages()[0].setPageSize(QgsLayoutSize( width, height, QgsUnitTypes.LayoutInches))
it exports without cropping but this time its width and height are not fixed as i specified


texts:
I have a WMS Layer published by a geoserver instance. It has the capability of being filtered by a CQL Filter via GetMap CQL_FILTER query parameter.
Adding it to QGIS works, but I can't find a place to add a CQL Filter in QGIS.
There is an old issue about it: https://github.com/qgis/QGIS/issues/15626
Is there any progress on the theme or a work around?
QGIS is version 3.10.4 - but I'm willing to change it if it solves my problem


texts: QGIS application is not displayed correctly & does not react

texts:
I'm having trouble with QGIS. I have tried many installations of it (removing the old one each time) and they are all having the same issue. It may be my computer but I don't even know what to look for. All of my other programs run fine.

I tried QGIS 2.14 32 bit, 2.14 64 bit, 2.16 32 bit, 2.16 64 bit, 2.18 32 bit, and 2.18 64 bit.

The problem is that when I go to click on any menus, drop downs, or buttons, either nothing happens or it shows only the edges of the drop down - no text or backgrounds. When I try to type in boxes, nothing happens. When I try to click a button to do something, nothing happens. The most I've ever been able to do was to add one vector layer with one field. Anything additional to that starts the problems.

I've rebooted my computer several times. QGIS is the only thing running. There are no other windows open. My computer is only ever using about 8% of its full potential at any one time running QGIS, so it's not "hanging".

I've tried my wired mouse and my touch pad. Both are doing the same thing. No problems with my mouse in any other program. Just verified all the drivers are up to date on both mice.

Any ideas of how to fix the issue?

Computer specs:

Windows 10 Home
Processor - Intel Core i7-6700HQ CPU @ 2.6 GHz
Graphics - NVIDIA GeForce® GTX 980M
16 GB RAM
System type - 64 GB


texts: PyQGIS grass7:v.distance parameter error "The number of columns and the number of upload parameters should be equal!"

texts:
I am currently working on a Python script for QGIS in which I would like to obtain the (minimum) distance from each point (polygon centroid) in the points layer centroid_layer to the line layer nile_layer.

I unfortunately still don't understand how the parameters need to be specified in order to store the distance in a column of the 'FROM' layer, here: centroid_layer.

Error message:

  Traceback (most recent call last):
    File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/code.py", line 90, in runcode
      exec(code, self.locals)
    File "", line 1, in
    File "", line 255, in
    File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/tools/general.py", line 106, in run
      return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
    File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/core/Processing.py", line 137, in runAlgorithm
      raise QgsProcessingException(msg)
  _core.QgsProcessingException: Unable to execute algorithm
  The number of columns and the number of upload parameters should be equal!

Code snippet:

# Import required packages
from PyQt5.QtGui import *
from PyQt5.QtCore import *
import os
import getpass
import requests
from qgis.core import *
from qgis.utils import iface
import qgis.utils
from qgis import core
import processing
from processing.core.Processing import Processing

[...]

# Distance to Nile
params = {
    'FROM' : centroid_layer , # 'From' layer, here: Murdock ethnicity centroid layer
    'FROM_TYPE' : 0 , # Type of 'from' layer; 0: point, 1: line, 2: boundary, 3: area, 4: centroid
    'TO': nile_layer,
    'TO_TYPE' : 1 , # Type of 'to' layer; 0: point, 1: line, 2: boundary, 3: area, 4: centroid
    'DMAX' : -1.0 , # Minimum distance or -1.0 for no limit
    'DMIN' : -1.0 , # Minimum distance or -1.0 for no limit
    'UPLOAD' : 1 , # 'dist' distance result type
    'COLUMN' : 'distance_to_nile' ,
    'OUTPUT' : 'Temporary datasets/Layers/distance_to_nile',
    '--OVERWRITE' : True
}
distance_nile_calculation = processing.run('grass7:v.distance', params)

I am running the code from the Python console within QGIS 3.10.0 on MacOS.


texts:
I have seen that OGR supports it's own form of SQL, I was wondering if it's possible to do spatial joins.

Something like select all objects where one intersects the other.  I did have a look around but couldn't really see anything.


texts: Geoserver - Access is denied to preview style

texts:
I am using Geoserver 2.10 and Tomcat 8.5.11. When I want edit style and preview it, geoserver give me following error:

  File failed to upload: Processing of multipart/form-data request failed. C:\Program Files\Apache Software Foundation\Tomcat 8.5\temp\upload__73bb119_15a50da8e4e__7fff_00000219.tmp (Access is denied)

Can you tell me how can I fix it ?

Here is image of error:


texts:
Is it possible to open local PDFs from within ESRI's AppStudio?  I can easily open pdfs and files from web links with "Qt.openUrlExternally('http..." but I need this app to work completely offline.  I've tried the approaches here with no luck:  Files and folders

Ultimately, I would want the AppStudio app to appear in the files app, so I could drop new or existing pdfs there via ITunes.  Is there a trick to getting the app to appear here, then how would I read/open files from that location?

Maybe there is a better approach?


texts: Why do I get a ". not nullable field error" when attempting to split a feature?

texts:
I'm trying to split feature, by point, into 2 new features.

 var pFeatureEdit = lineToSplit.Feature as IFeatureEdit;
 var set = pFeatureEdit.Split(point);

When this feature has configured relation, I got exception.

Got next COM error: The field is not nullable [FIELDID]
Where FIELDID - not nullable foreign key in database

Anybody know how to solve this issue?

UPD. When I split this feature using standard Split Tool in ArcMap - it splits successfully. This uses another method to split features?


texts:
I added a MySQL layer containing GeoData to my QGIS 3 project.
I want this layer updated every time the database updates.

Is there a way to do that?


texts: Select by attributes with user input using Python

texts:
I am new to GIS and Python is not my forte, so I am looking for some help. Ideally, what I am trying to do is create a tool that allows the user to select an address from my Address points layer. I need to select values from multiple fields: STREET_NUMBER, STREET_NAME and STREET_TYPE. What I really need help with is building my expression. Here is what I have so far:

# Import system modules
import arcpy
from arcpy import env

# Set workspace
workspace = r"O:\Users\Student\Scratch\SC_Test.gdb"

# Set variables
stNum = arcpy.GetParameterAsText(0)
stName = arcpy.GetParameterAsText(1)
stType = arcpy.GetParameterAsText(2)

# Make layer from feature class
arcpy.MakeFeatureLayer_management("Address", "add_lyr")

arcpy.SelectLayerByAttribute_management("add_lyr", "NEW_SELECTION", "STREET_NUMBER" =  stNum)


texts: Exporting GeoPackage layer to GeoJSON changes CRS from WGS 84 EPSG:4326 to OGC:CRS84

texts:
I'm attempting to convert a polygon layer from GeoPackage format to GeoJSON in QGIS 3.10.1. The original layer has a coordinate reference system of WGS 84 (EPSG:4326) which is what I need. Yet the resulting GeoJSON output has a coordinate reference system of OGC:CRS84 which is causing invalid GeoJSON errors when attempting to import it into various web applications.

i.e. GeoJSONLint, Right Hand Rule GeoJSON Fixer, etc. I keep getting the "Line 1: old-style crs member is not recommended, this object is equivalent to the default and should be removed" error.

How can I get the resulting GeoJSON file to stop converting the CRS and keep WGS 84 (EPSG:4326). Or is this the intended behavior of QGIS upon export to GeoJSON and I'm thinking about this wrong?


texts: Working with address data extracted from *.osm

texts:
My goal is get all the addresses from a certain country, ideally - from all countries, but for now - only from one.

I downloaded an *.osm file from here http://download.geofabrik.de/ for a certain territory which contained a few countries and extracted the data from that file using Osm2pgsql to a Postgres database. There were no errors, all went well.

But now I can't get along with those tables in the db, I can't figure out how actually to query all the addresses from a country I want?


texts: get active colourramp form qgis rasterlayer

texts:
I m building a pythonscript were i need the active singleband pseudocolour colourband of a raster. I want to give the active colourramp a new min and max value.

I can only find examples like this one. In this example they create a new colourramp. But i want to keep the exsiting one.

Does somebody know how i can get the active colour ramp with python?


texts: OpenLayers 3 ol.source.xyz set custom origin

texts:
I am using local tiles from my file system in OpenLayers 3 (ol.source.xyz), and my local tiles are stored as  0/0/0.png and so on. Now the thing is that when I use these tiles on my map it starts from the top left corner (0,0), instead of I want to set custom origin of my tiles, is there any way to change the origin (instead of 0,0, I want to use my own custom coordinates).

Here is my code:


texts: Use 2 different symbologies on the same feature QGIS

texts:
I have data symbolized with different categories
Each category has a field for different symbology style like the table below.
Is it possible to visualize  the same category of data with in different styles without adding the layer twice?

I would like to do something like this
Category=Forest AND Symbology=1 THEN 'Green'
Category=Forest AND Symbology=2 THEN 'Blue'

In QGIS there is a selection "Data defined properties"  with an expression builder, could it be done from it maybe?

Category Symbology
Forest    1
Forest  2
Shrub   1
Shrub   2
Grass   1
Grass   2


texts: Zonal Histogram Output Problems - No values associated with the habitat type

texts:
I have a outline of some islands (shapefile or raster) and GAP landcover data (raster). I am trying to use Zonal Histogram to create a file with a histogram of the habitat types in each island. Each habitat type has a value that is the field I'm computing the histogram for.

Sounds easy enough, right? Well, ArcGIS 10.3 has problems with the zonal histogram, so I'm also using ArcGIS 9.3 SP1 on a virtual machine.

So here's the problem - the file that ArcGIS creates doesn't include that habitat ID value in the histogram - the rows are just labeled 1-255 sequentially  (I have 115 habitat types). I can't figure out what order they are actually in and I can't make any sense of how to paste them in by hand. There are a number of lines that are just 0's across the board. Deleting the 0 rows does not result in 115 rows.

I've been plugging away at this in both versions of ArcGIS and in R for about two weeks now and keep hitting this wall.


texts: How can a style be set on a specific polygon using OpenLayers?

texts:
Using a polygon id, I want to highlight a polygon, while the mouse is outside of the map area. I can identify the polygon using the following approach:

x = 53;
for(i=0; i<vector_layer.features.length; i++){
        y = vector_layer.features[i].attributes.id
        if(x == y){
            vector_layer.features[i].style = ??
        }
}

However, I can't figure out how to highlight the polygon even though I know which one it is. I have previously defined a style for when the mouse hovers over the polygon and it is selected, but I don't know how to apply this style to the chosen polygon.


texts: Exporting layers to different sheets of an Excel file

texts:
I have two memory layers that I would like to write their attribute tables to the same Excel file but in different sheets. So far, I've managed to write the attribute tables to separate files using QgsVectorFileWriter.writeAsVectorFormat.
layer1 = QgsProject.instance().mapLayersByName("memory_layer1")[0]
layer2 = QgsProject.instance().mapLayersByName("memory_layer2")[0]

# export to separate file
QgsVectorFileWriter.writeAsVectorFormat(layer1, 'path/file1.xlsx',
    "utf-8", layer1.crs(), "XLSX", layerOptions = ['GEOMETRY=AS_XYZ'])

QgsVectorFileWriter.writeAsVectorFormat(layer2, 'path/file2.xlsx',
    "utf-8", layer2.crs(), "XLSX", layerOptions = ['GEOMETRY=AS_XYZ'])


texts: Spatial join in R - Adding points to polygons with multiple corresponding points

texts:
I have carried out a spatial join in R in order to join a points shapefile to the corresponding polygon in a polygon shapefile using over() from the sp package and spCbind().

join <- over(polygons, points)

join <- spCbind(polygons, join)

However, some of the polygons have multiple corresponding points. Using this method, only one of these corresponding points is maintained in the joined dataset.

Is there a method for carrying out this task that will keep the data of all of the corresponding points for each polygon in the resulting dataset?

I would prefer an R solution if possible but I can also use QGIS or ArcGIS if they have tools to solve this.


texts: Doing weighted overlay analysis using CSV data joined with shapefile in QGIS?

texts:
I have a CSV file which has data for different counties income, Population and Proximity of United States of America like this :

]1

I have already joined the CSV file with the shapefile of United States map and duplicated them to use them as different layers.

How can I use them to do weighted overlay analysis to find the list of top 10 counties with respect to these 3 attributes and respective weights?


texts: Moving QGIS basemap by hand in order to achieve overlaying with point layer?

texts:
I am an environmental engineer and the aim of my QGIS 2.0 project is to visualize the results of a  sewage water system inspection.
In a first step I included the hoppers from a csv file and added a google streets layer.

As you can see in the picture the hoppers don't overlay with the base map. I compared to the coordinate calculation of a mountain climber page (http://www.deine-berge.de/Rechner/Koordinaten/Sandersdorf,-06792-Sandersdorf-Brehna,-Deutschland), which confirms the on the fly crs transformation of qgis 2.0.
So the coordinates of the hoppers are may wrong. However it is the only csv data I have and in the GIS system of our branch office it fits to the base map.
Is there a possibility to move the base map by hand in order to achieve an overlaying?
For my purpose I don't need an exact fitting, but a difference of 500 m is too much.


texts: Outputting same name as input using Copy Features in ModelBuilder?

texts:
I am using ModelBuilder to iterate Closest Facility analysis for every building at an archaeological site (OD Cost Matrix is not what I need--I need the geometry output, not just a cost matrix table). My draft model is below. One question: at the end of the model, it takes the routes generated from the Closest Facility analysis and copies it to a new feature class. However, the Copy Features tool seems to require a specific name for the output. When I iterate to the next facility, it will either overwrite or break the model.

So the question is simple: how do I configure Copy Features tool to output the name from the corresponding iteration (taken form the ObjectID field from the Iterate Feature Selection)?


texts: Can't find function st_join in R sf package

texts:
I am trying to run the st_join example from the sf user manual https://cran.r-project.org/web/packages/sf/sf.pdf (p. 43).

a = st_sf(a = 1:3, geom = st_sfc(st_point(c(1,1)), st_point(c(2,2)), st_point(c(3,3))))

b = st_sf(a = 11:14, geom = st_sfc(st_point(c(10,10)), st_point(c(2,2)), st_point(c(2,2)), st_point(c(3,3))))

st_join(a, b)

I get the error message that R can't find the function sf_join.
I am running R 3.3.2 on x86_64-w64-mingw32
I installed sf0.3.1

What am I doing wrong?